diff --git a/.gitignore b/.gitignore index 723e9f7..ceb1e40 100644 --- a/.gitignore +++ b/.gitignore @@ -4,4 +4,6 @@ *tar.gz +*out + TEST/* \ No newline at end of file diff --git a/ABACUS.org b/ABACUS.org index cf3b5f0..2cb6bdd 100644 --- a/ABACUS.org +++ b/ABACUS.org @@ -32,15 +32,6 @@ Type your description here TIP: Search for the string BUGRISK in the codebase -** BUGRISK Value of LiebLin ln_Density_ME - - State "BUGRISK" from "" [2018-02-11 Sun 09:11] - -File: ln_Density_ME.cc -line 66 - -Why real? - - * Priority :ABACUS:Dev:Priority: :PROPERTIES: :ARCHIVE: %s_archive::* Priority diff --git a/ABACUS.org_archive b/ABACUS.org_archive new file mode 100644 index 0000000..0eaece9 --- /dev/null +++ b/ABACUS.org_archive @@ -0,0 +1,26 @@ +# -*- mode: org -*- + + +Archived entries from file /Users/jscaux/WORK/_ABACUS/ABACUS/ABACUS.org + + +* Bugrisks + +** SOLVED Value of LiebLin ln_Density_ME + CLOSED: [2018-02-15 Thu 08:42] + - State "SOLVED" from "BUGRISK" [2018-02-15 Thu 08:42] \\ + No problem. Product is always real. + - State "BUGRISK" from "" [2018-02-11 Sun 09:11] + :PROPERTIES: + :ARCHIVE_TIME: 2018-02-15 Thu 08:42 + :ARCHIVE_FILE: ~/WORK/_ABACUS/ABACUS/ABACUS.org + :ARCHIVE_OLPATH: Bugrisks + :ARCHIVE_CATEGORY: ABACUS + :ARCHIVE_TODO: SOLVED + :ARCHIVE_ITAGS: ABACUS Dev Bugrisks + :END: + +File: ln_Density_ME.cc +line 66 + +Why real? diff --git a/ABACUS_Doxyfile b/ABACUS_Doxyfile index 62de073..558bf0a 100644 --- a/ABACUS_Doxyfile +++ b/ABACUS_Doxyfile @@ -38,7 +38,7 @@ PROJECT_NAME = ABACUS # could be handy for archiving the generated documentation or if some version # control system is used. -PROJECT_NUMBER = +PROJECT_NUMBER = # Using the PROJECT_BRIEF tag one can provide an optional one line description # for a project that appears at the top of each page and should give viewer a @@ -51,7 +51,7 @@ PROJECT_BRIEF = "Numerical tools for integrable models" # pixels and the maximum width should not exceed 200 pixels. Doxygen will copy # the logo to the output directory. -PROJECT_LOGO = +PROJECT_LOGO = # The OUTPUT_DIRECTORY tag is used to specify the (relative or absolute) path # into which the generated documentation will be written. If a relative path is @@ -162,7 +162,7 @@ FULL_PATH_NAMES = YES # will be relative from the directory where doxygen is started. # This tag requires that the tag FULL_PATH_NAMES is set to YES. -STRIP_FROM_PATH = +STRIP_FROM_PATH = # The STRIP_FROM_INC_PATH tag can be used to strip a user-defined part of the # path mentioned in the documentation of a class, which tells the reader which @@ -171,7 +171,7 @@ STRIP_FROM_PATH = # specify the list of include paths that are normally passed to the compiler # using the -I flag. -STRIP_FROM_INC_PATH = +STRIP_FROM_INC_PATH = # If the SHORT_NAMES tag is set to YES, doxygen will generate much shorter (but # less readable) file names. This can be useful is your file systems doesn't @@ -238,13 +238,13 @@ TAB_SIZE = 4 # "Side Effects:". You can put \n's in the value part of an alias to insert # newlines. -ALIASES = +ALIASES = # This tag can be used to specify a number of word-keyword mappings (TCL only). # A mapping has the form "name=value". For example adding "class=itcl::class" # will allow you to use the command class in the itcl::class meaning. -TCL_SUBST = +TCL_SUBST = # Set the OPTIMIZE_OUTPUT_FOR_C tag to YES if your project consists of C sources # only. Doxygen will then generate output that is more tailored for C. For @@ -291,7 +291,7 @@ OPTIMIZE_OUTPUT_VHDL = NO # Note that for custom extensions you also need to set FILE_PATTERNS otherwise # the files are not read by doxygen. -EXTENSION_MAPPING = +EXTENSION_MAPPING = # If the MARKDOWN_SUPPORT tag is enabled then doxygen pre-processes all comments # according to the Markdown format, which allows for more readable @@ -648,7 +648,7 @@ GENERATE_DEPRECATEDLIST= YES # sections, marked by \if ... \endif and \cond # ... \endcond blocks. -ENABLED_SECTIONS = +ENABLED_SECTIONS = # The MAX_INITIALIZER_LINES tag determines the maximum number of lines that the # initial value of a variable or macro / define can have for it to appear in the @@ -690,7 +690,7 @@ SHOW_NAMESPACES = YES # by doxygen. Whatever the program writes to standard output is used as the file # version. For an example see the documentation. -FILE_VERSION_FILTER = +FILE_VERSION_FILTER = # The LAYOUT_FILE tag can be used to specify a layout file which will be parsed # by doxygen. The layout file controls the global structure of the generated @@ -703,7 +703,7 @@ FILE_VERSION_FILTER = # DoxygenLayout.xml, doxygen will parse it automatically even if the LAYOUT_FILE # tag is left empty. -LAYOUT_FILE = +LAYOUT_FILE = # The CITE_BIB_FILES tag can be used to specify one or more bib files containing # the reference definitions. This must be a list of .bib files. The .bib @@ -713,7 +713,7 @@ LAYOUT_FILE = # LATEX_BIB_STYLE. To use this feature you need bibtex and perl available in the # search path. See also \cite for info how to create references. -CITE_BIB_FILES = +CITE_BIB_FILES = #--------------------------------------------------------------------------- # Configuration options related to warning and progress messages @@ -778,7 +778,7 @@ WARN_FORMAT = "$file:$line: $text" # messages should be written. If left blank the output is written to standard # error (stderr). -WARN_LOGFILE = +WARN_LOGFILE = #--------------------------------------------------------------------------- # Configuration options related to the input files @@ -873,7 +873,7 @@ RECURSIVE = YES # Note that relative paths are relative to the directory from which doxygen is # run. -EXCLUDE = +EXCLUDE = # The EXCLUDE_SYMLINKS tag can be used to select whether or not files or # directories that are symbolic links (a Unix file system feature) are excluded @@ -889,7 +889,7 @@ EXCLUDE_SYMLINKS = NO # Note that the wildcards are matched against the file with absolute path, so to # exclude all test directories for example use the pattern */test/* -EXCLUDE_PATTERNS = +EXCLUDE_PATTERNS = # The EXCLUDE_SYMBOLS tag can be used to specify one or more symbol names # (namespaces, classes, functions, etc.) that should be excluded from the @@ -900,13 +900,13 @@ EXCLUDE_PATTERNS = # Note that the wildcards are matched against the file with absolute path, so to # exclude all test directories use the pattern */test/* -EXCLUDE_SYMBOLS = +EXCLUDE_SYMBOLS = # The EXAMPLE_PATH tag can be used to specify one or more files or directories # that contain example code fragments that are included (see the \include # command). -EXAMPLE_PATH = +EXAMPLE_PATH = # If the value of the EXAMPLE_PATH tag contains directories, you can use the # EXAMPLE_PATTERNS tag to specify one or more wildcard pattern (like *.cpp and @@ -926,7 +926,7 @@ EXAMPLE_RECURSIVE = NO # that contain images that are to be included in the documentation (see the # \image command). -IMAGE_PATH = +IMAGE_PATH = # The INPUT_FILTER tag can be used to specify a program that doxygen should # invoke to filter for each input file. Doxygen will invoke the filter program @@ -947,7 +947,7 @@ IMAGE_PATH = # need to set EXTENSION_MAPPING for the extension otherwise the files are not # properly processed by doxygen. -INPUT_FILTER = +INPUT_FILTER = # The FILTER_PATTERNS tag can be used to specify filters on a per file pattern # basis. Doxygen will compare the file name with each pattern and apply the @@ -960,7 +960,7 @@ INPUT_FILTER = # need to set EXTENSION_MAPPING for the extension otherwise the files are not # properly processed by doxygen. -FILTER_PATTERNS = +FILTER_PATTERNS = # If the FILTER_SOURCE_FILES tag is set to YES, the input filter (if set using # INPUT_FILTER) will also be used to filter the input files that are used for @@ -975,14 +975,14 @@ FILTER_SOURCE_FILES = NO # *.ext= (so without naming a filter). # This tag requires that the tag FILTER_SOURCE_FILES is set to YES. -FILTER_SOURCE_PATTERNS = +FILTER_SOURCE_PATTERNS = # If the USE_MDFILE_AS_MAINPAGE tag refers to the name of a markdown file that # is part of the input, its contents will be placed on the main page # (index.html). This can be useful if you have a project on for instance GitHub # and want to reuse the introduction page also for the doxygen output. -USE_MDFILE_AS_MAINPAGE = +USE_MDFILE_AS_MAINPAGE = #--------------------------------------------------------------------------- # Configuration options related to source browsing @@ -1087,7 +1087,7 @@ CLANG_ASSISTED_PARSING = NO # specified with INPUT and INCLUDE_PATH. # This tag requires that the tag CLANG_ASSISTED_PARSING is set to YES. -CLANG_OPTIONS = +CLANG_OPTIONS = #--------------------------------------------------------------------------- # Configuration options related to the alphabetical class index @@ -1113,7 +1113,7 @@ COLS_IN_ALPHA_INDEX = 5 # while generating the index headers. # This tag requires that the tag ALPHABETICAL_INDEX is set to YES. -IGNORE_PREFIX = +IGNORE_PREFIX = #--------------------------------------------------------------------------- # Configuration options related to the HTML output @@ -1157,7 +1157,7 @@ HTML_FILE_EXTENSION = .html # of the possible markers and block names see the documentation. # This tag requires that the tag GENERATE_HTML is set to YES. -HTML_HEADER = +HTML_HEADER = # The HTML_FOOTER tag can be used to specify a user-defined HTML footer for each # generated HTML page. If the tag is left blank doxygen will generate a standard @@ -1167,7 +1167,7 @@ HTML_HEADER = # that doxygen normally uses. # This tag requires that the tag GENERATE_HTML is set to YES. -HTML_FOOTER = +HTML_FOOTER = # The HTML_STYLESHEET tag can be used to specify a user-defined cascading style # sheet that is used by each HTML page. It can be used to fine-tune the look of @@ -1179,7 +1179,7 @@ HTML_FOOTER = # obsolete. # This tag requires that the tag GENERATE_HTML is set to YES. -HTML_STYLESHEET = +HTML_STYLESHEET = # The HTML_EXTRA_STYLESHEET tag can be used to specify additional user-defined # cascading style sheets that are included after the standard style sheets @@ -1192,7 +1192,7 @@ HTML_STYLESHEET = # list). For an example see the documentation. # This tag requires that the tag GENERATE_HTML is set to YES. -HTML_EXTRA_STYLESHEET = +HTML_EXTRA_STYLESHEET = # The HTML_EXTRA_FILES tag can be used to specify one or more extra images or # other source files which should be copied to the HTML output directory. Note @@ -1202,7 +1202,7 @@ HTML_EXTRA_STYLESHEET = # files will be copied as-is; there are no commands or markers available. # This tag requires that the tag GENERATE_HTML is set to YES. -HTML_EXTRA_FILES = +HTML_EXTRA_FILES = # The HTML_COLORSTYLE_HUE tag controls the color of the HTML output. Doxygen # will adjust the colors in the style sheet and background images according to @@ -1331,7 +1331,7 @@ GENERATE_HTMLHELP = NO # written to the html output directory. # This tag requires that the tag GENERATE_HTMLHELP is set to YES. -CHM_FILE = +CHM_FILE = # The HHC_LOCATION tag can be used to specify the location (absolute path # including file name) of the HTML help compiler (hhc.exe). If non-empty, @@ -1339,7 +1339,7 @@ CHM_FILE = # The file has to be specified with full path. # This tag requires that the tag GENERATE_HTMLHELP is set to YES. -HHC_LOCATION = +HHC_LOCATION = # The GENERATE_CHI flag controls if a separate .chi index file is generated # (YES) or that it should be included in the master .chm file (NO). @@ -1352,7 +1352,7 @@ GENERATE_CHI = NO # and project file content. # This tag requires that the tag GENERATE_HTMLHELP is set to YES. -CHM_INDEX_ENCODING = +CHM_INDEX_ENCODING = # The BINARY_TOC flag controls whether a binary table of contents is generated # (YES) or a normal table of contents (NO) in the .chm file. Furthermore it @@ -1383,7 +1383,7 @@ GENERATE_QHP = NO # the HTML output folder. # This tag requires that the tag GENERATE_QHP is set to YES. -QCH_FILE = +QCH_FILE = # The QHP_NAMESPACE tag specifies the namespace to use when generating Qt Help # Project output. For more information please see Qt Help Project / Namespace @@ -1408,7 +1408,7 @@ QHP_VIRTUAL_FOLDER = doc # filters). # This tag requires that the tag GENERATE_QHP is set to YES. -QHP_CUST_FILTER_NAME = +QHP_CUST_FILTER_NAME = # The QHP_CUST_FILTER_ATTRS tag specifies the list of the attributes of the # custom filter to add. For more information please see Qt Help Project / Custom @@ -1416,21 +1416,21 @@ QHP_CUST_FILTER_NAME = # filters). # This tag requires that the tag GENERATE_QHP is set to YES. -QHP_CUST_FILTER_ATTRS = +QHP_CUST_FILTER_ATTRS = # The QHP_SECT_FILTER_ATTRS tag specifies the list of the attributes this # project's filter section matches. Qt Help Project / Filter Attributes (see: # http://qt-project.org/doc/qt-4.8/qthelpproject.html#filter-attributes). # This tag requires that the tag GENERATE_QHP is set to YES. -QHP_SECT_FILTER_ATTRS = +QHP_SECT_FILTER_ATTRS = # The QHG_LOCATION tag can be used to specify the location of Qt's # qhelpgenerator. If non-empty doxygen will try to run qhelpgenerator on the # generated .qhp file. # This tag requires that the tag GENERATE_QHP is set to YES. -QHG_LOCATION = +QHG_LOCATION = # If the GENERATE_ECLIPSEHELP tag is set to YES, additional index files will be # generated, together with the HTML files, they form an Eclipse help plugin. To @@ -1533,7 +1533,7 @@ FORMULA_TRANSPARENT = YES # The default value is: NO. # This tag requires that the tag GENERATE_HTML is set to YES. -USE_MATHJAX = NO +USE_MATHJAX = YES # When MathJax is enabled you can set the default output format to be used for # the MathJax output. See the MathJax site (see: @@ -1563,7 +1563,7 @@ MATHJAX_RELPATH = http://cdn.mathjax.org/mathjax/latest # MATHJAX_EXTENSIONS = TeX/AMSmath TeX/AMSsymbols # This tag requires that the tag USE_MATHJAX is set to YES. -MATHJAX_EXTENSIONS = +MATHJAX_EXTENSIONS = # The MATHJAX_CODEFILE tag can be used to specify a file with javascript pieces # of code that will be used on startup of the MathJax code. See the MathJax site @@ -1571,7 +1571,7 @@ MATHJAX_EXTENSIONS = # example see the documentation. # This tag requires that the tag USE_MATHJAX is set to YES. -MATHJAX_CODEFILE = +MATHJAX_CODEFILE = # When the SEARCHENGINE tag is enabled doxygen will generate a search box for # the HTML output. The underlying search engine uses javascript and DHTML and @@ -1631,7 +1631,7 @@ EXTERNAL_SEARCH = NO # Searching" for details. # This tag requires that the tag SEARCHENGINE is set to YES. -SEARCHENGINE_URL = +SEARCHENGINE_URL = # When SERVER_BASED_SEARCH and EXTERNAL_SEARCH are both enabled the unindexed # search data is written to a file for indexing by an external tool. With the @@ -1647,7 +1647,7 @@ SEARCHDATA_FILE = searchdata.xml # projects and redirect the results back to the right project. # This tag requires that the tag SEARCHENGINE is set to YES. -EXTERNAL_SEARCH_ID = +EXTERNAL_SEARCH_ID = # The EXTRA_SEARCH_MAPPINGS tag can be used to enable searching through doxygen # projects other than the one defined by this configuration file, but that are @@ -1657,7 +1657,7 @@ EXTERNAL_SEARCH_ID = # EXTRA_SEARCH_MAPPINGS = tagname1=loc1 tagname2=loc2 ... # This tag requires that the tag SEARCHENGINE is set to YES. -EXTRA_SEARCH_MAPPINGS = +EXTRA_SEARCH_MAPPINGS = #--------------------------------------------------------------------------- # Configuration options related to the LaTeX output @@ -1721,7 +1721,7 @@ PAPER_TYPE = a4 # If left blank no extra packages will be included. # This tag requires that the tag GENERATE_LATEX is set to YES. -EXTRA_PACKAGES = +EXTRA_PACKAGES = # The LATEX_HEADER tag can be used to specify a personal LaTeX header for the # generated LaTeX document. The header should contain everything until the first @@ -1737,7 +1737,7 @@ EXTRA_PACKAGES = # to HTML_HEADER. # This tag requires that the tag GENERATE_LATEX is set to YES. -LATEX_HEADER = +LATEX_HEADER = # The LATEX_FOOTER tag can be used to specify a personal LaTeX footer for the # generated LaTeX document. The footer should contain everything after the last @@ -1748,7 +1748,7 @@ LATEX_HEADER = # Note: Only use a user-defined footer if you know what you are doing! # This tag requires that the tag GENERATE_LATEX is set to YES. -LATEX_FOOTER = +LATEX_FOOTER = # The LATEX_EXTRA_STYLESHEET tag can be used to specify additional user-defined # LaTeX style sheets that are included after the standard style sheets created @@ -1759,7 +1759,7 @@ LATEX_FOOTER = # list). # This tag requires that the tag GENERATE_LATEX is set to YES. -LATEX_EXTRA_STYLESHEET = +LATEX_EXTRA_STYLESHEET = # The LATEX_EXTRA_FILES tag can be used to specify one or more extra images or # other source files which should be copied to the LATEX_OUTPUT output @@ -1767,7 +1767,7 @@ LATEX_EXTRA_STYLESHEET = # markers available. # This tag requires that the tag GENERATE_LATEX is set to YES. -LATEX_EXTRA_FILES = +LATEX_EXTRA_FILES = # If the PDF_HYPERLINKS tag is set to YES, the LaTeX that is generated is # prepared for conversion to PDF (using ps2pdf or pdflatex). The PDF file will @@ -1875,14 +1875,14 @@ RTF_HYPERLINKS = NO # default style sheet that doxygen normally uses. # This tag requires that the tag GENERATE_RTF is set to YES. -RTF_STYLESHEET_FILE = +RTF_STYLESHEET_FILE = # Set optional variables used in the generation of an RTF document. Syntax is # similar to doxygen's config file. A template extensions file can be generated # using doxygen -e rtf extensionFile. # This tag requires that the tag GENERATE_RTF is set to YES. -RTF_EXTENSIONS_FILE = +RTF_EXTENSIONS_FILE = # If the RTF_SOURCE_CODE tag is set to YES then doxygen will include source code # with syntax highlighting in the RTF output. @@ -1927,7 +1927,7 @@ MAN_EXTENSION = .3 # MAN_EXTENSION with the initial . removed. # This tag requires that the tag GENERATE_MAN is set to YES. -MAN_SUBDIR = +MAN_SUBDIR = # If the MAN_LINKS tag is set to YES and doxygen generates man output, then it # will generate one additional man file for each entity documented in the real @@ -2040,7 +2040,7 @@ PERLMOD_PRETTY = YES # overwrite each other's variables. # This tag requires that the tag GENERATE_PERLMOD is set to YES. -PERLMOD_MAKEVAR_PREFIX = +PERLMOD_MAKEVAR_PREFIX = #--------------------------------------------------------------------------- # Configuration options related to the preprocessor @@ -2081,7 +2081,7 @@ SEARCH_INCLUDES = YES # preprocessor. # This tag requires that the tag SEARCH_INCLUDES is set to YES. -INCLUDE_PATH = +INCLUDE_PATH = # You can use the INCLUDE_FILE_PATTERNS tag to specify one or more wildcard # patterns (like *.h and *.hpp) to filter out the header-files in the @@ -2089,7 +2089,7 @@ INCLUDE_PATH = # used. # This tag requires that the tag ENABLE_PREPROCESSING is set to YES. -INCLUDE_FILE_PATTERNS = +INCLUDE_FILE_PATTERNS = # The PREDEFINED tag can be used to specify one or more macro names that are # defined before the preprocessor is started (similar to the -D option of e.g. @@ -2099,7 +2099,7 @@ INCLUDE_FILE_PATTERNS = # recursively expanded use the := operator instead of the = operator. # This tag requires that the tag ENABLE_PREPROCESSING is set to YES. -PREDEFINED = +PREDEFINED = # If the MACRO_EXPANSION and EXPAND_ONLY_PREDEF tags are set to YES then this # tag can be used to specify a list of macro names that should be expanded. The @@ -2108,7 +2108,7 @@ PREDEFINED = # definition found in the source code. # This tag requires that the tag ENABLE_PREPROCESSING is set to YES. -EXPAND_AS_DEFINED = +EXPAND_AS_DEFINED = # If the SKIP_FUNCTION_MACROS tag is set to YES then doxygen's preprocessor will # remove all references to function-like macros that are alone on a line, have @@ -2137,13 +2137,13 @@ SKIP_FUNCTION_MACROS = YES # the path). If a tag file is not located in the directory in which doxygen is # run, you must also specify the path to the tagfile here. -TAGFILES = +TAGFILES = # When a file name is specified after GENERATE_TAGFILE, doxygen will create a # tag file that is based on the input files it reads. See section "Linking to # external documentation" for more information about the usage of tag files. -GENERATE_TAGFILE = +GENERATE_TAGFILE = # If the ALLEXTERNALS tag is set to YES, all external class will be listed in # the class index. If set to NO, only the inherited external classes will be @@ -2192,14 +2192,14 @@ CLASS_DIAGRAMS = YES # the mscgen tool resides. If left empty the tool is assumed to be found in the # default search path. -MSCGEN_PATH = +MSCGEN_PATH = # You can include diagrams made with dia in doxygen documentation. Doxygen will # then run dia to produce the diagram and insert it in the documentation. The # DIA_PATH tag allows you to specify the directory where the dia binary resides. # If left empty dia is assumed to be found in the default search path. -DIA_PATH = +DIA_PATH = # If set to YES the inheritance and collaboration graphs will hide inheritance # and usage relations if the target is undocumented or is not a class. @@ -2248,7 +2248,7 @@ DOT_FONTSIZE = 10 # the path where dot can find it using this tag. # This tag requires that the tag HAVE_DOT is set to YES. -DOT_FONTPATH = +DOT_FONTPATH = # If the CLASS_GRAPH tag is set to YES then doxygen will generate a graph for # each documented class showing the direct and indirect inheritance relations. @@ -2392,26 +2392,26 @@ INTERACTIVE_SVG = NO # found. If left blank, it is assumed the dot tool can be found in the path. # This tag requires that the tag HAVE_DOT is set to YES. -DOT_PATH = +DOT_PATH = # The DOTFILE_DIRS tag can be used to specify one or more directories that # contain dot files that are included in the documentation (see the \dotfile # command). # This tag requires that the tag HAVE_DOT is set to YES. -DOTFILE_DIRS = +DOTFILE_DIRS = # The MSCFILE_DIRS tag can be used to specify one or more directories that # contain msc files that are included in the documentation (see the \mscfile # command). -MSCFILE_DIRS = +MSCFILE_DIRS = # The DIAFILE_DIRS tag can be used to specify one or more directories that # contain dia files that are included in the documentation (see the \diafile # command). -DIAFILE_DIRS = +DIAFILE_DIRS = # When using plantuml, the PLANTUML_JAR_PATH tag should be used to specify the # path where java can find the plantuml.jar file. If left blank, it is assumed @@ -2419,17 +2419,17 @@ DIAFILE_DIRS = # generate a warning when it encounters a \startuml command in this case and # will not generate output for the diagram. -PLANTUML_JAR_PATH = +PLANTUML_JAR_PATH = # When using plantuml, the PLANTUML_CFG_FILE tag can be used to specify a # configuration file for plantuml. -PLANTUML_CFG_FILE = +PLANTUML_CFG_FILE = # When using plantuml, the specified paths are searched for files specified by # the !include statement in a plantuml block. -PLANTUML_INCLUDE_PATH = +PLANTUML_INCLUDE_PATH = # The DOT_GRAPH_MAX_NODES tag can be used to set the maximum number of nodes # that will be shown in the graph. If the number of nodes in a graph becomes diff --git a/ABACUS_Usage_Example_Heis.cc b/ABACUS_Usage_Example_Heis.cc new file mode 100644 index 0000000..ef33418 --- /dev/null +++ b/ABACUS_Usage_Example_Heis.cc @@ -0,0 +1,125 @@ +/********************************************************** + +This software is part of J.-S. Caux's ABACUS++ library. + +Copyright (c) J.-S. Caux + +----------------------------------------------------------- + +File: ABACUS+_Usage_Example_Heis.cc + +Purpose: illustrates basic use of ABACUS for spin chains. + +***********************************************************/ + +#include "ABACUS.h" + +using namespace std; +using namespace JSC; + + +int main() +{ + clock_t StartTime = clock(); + + // Refer to include/JSC_Heis.h for all class definitions + // and to src/HEIS/ for the actual implementations. + + // Set basic system parameters: + DP Delta = 1.0; // anisotropy + int N = 64; // chain length + int M = 32; // number of downturned spins as compared to all spins up; must be <= N/2 (on or above equator) + + // Define the chain: J, Delta, h, Nsites + Heis_Chain chain(1.0, Delta, 0.0, N); + + // The Heis_Chain class so constructed contains information about all types of strings: + cout << "Delta = " << Delta << endl << "Possible string lengths and parities: "; + for (int j = 0; j < chain.Nstrings; ++j) cout << "(" << chain.Str_L[j] << ", " << chain.par[j] << ")" << "\t"; + cout << endl; + + + // Constructing a Bethe state: start with the ground state for simplicity. + // Define the base: chain, Mdown + // A base by definition represents a partition of the M down spins into different strings. + Heis_Base gbase(chain, M); // This puts all the M down spins into one-strings. + + // Define the ground state + //XXZ_Bethe_State gstate(chain, gbase); // Use this constructor for 0 < Delta < 1 + XXX_Bethe_State gstate(chain, gbase); // Use this constructor for Delta == 1 + //XXZ_gpd_Bethe_State gstate(chain, gbase); // Use this constructor for Delta > 1. + // Anisotropies Delta < 0 are not separately implemented. + + // The assignment operator is overloaded for Bethe states: + XXX_Bethe_State gstate2 = gstate; // copies all data of gstate into new object gstate2 + + + // Compute everything about the ground state + gstate.Compute_All(true); + + // Output information about the state + cout << gstate << endl; + + + // Now define some excited state. + + // First method: + // Start with defining a base using a base constructor with a vector of numbers of rapidities of each type as argument: + // First define Nrapidities vector: + Vect Nrapidities(0, chain.Nstrings); + // Assigning the numbers that follow, you have to ensure (yourself!) the \sum_j M_j n_j = M (where n_j is the length of type j). + Nrapidities[0] = M-2; // number of one-strings + Nrapidities[1] = 1; // one two-string + // Define the base: + Heis_Base ebase(chain, Nrapidities); + // Once the base is defined, the limiting quantum numbers are automatically computed: + cout << "ebase defined, data is (Mdown, Nrap, Nraptot, Ix2_infty, Ix2_min, Ix2_max, baselabel):" + << ebase.Mdown << endl << ebase.Nrap << endl << ebase.Nraptot << endl << ebase.Ix2_infty << endl + << ebase.Ix2_min << endl << ebase.Ix2_max << endl << ebase.baselabel << endl; + + // An excited state can then be defined using this new base: + XXX_Bethe_State estate(chain, ebase); + // Individual quantum numbers can be manipulated: this will NOT update the state label or verify range of Ix2 + estate.Ix2[0][0] = M+1; + estate.Compute_All(true); + cout << endl << "estate: " << estate << endl; + + + // Second method of constructing a base: + Heis_Base ebase2(chain, "31"); // This is another constructor for all rapidities in one-strings + //Heis_Base ebase(chain, "29x1y1"); // one two-string (XXX) + + + // Construct a new Bethe state + XXX_Bethe_State estateref (chain, ebase2); // this will contain the lowest-energy quantum numbers for this base + XXX_Bethe_State estate2 (chain, ebase2); // yet another state + + // Setting a state to a given label: + estate2.Set_to_Label ("31_1_nh", estateref.Ix2); // The base of estate must coincide with the base in the label. Label is relative to estateref.Ix2 + estate2.Compute_All(true); + cout << "estate2: " << estate2 << endl; + + // Energy and momentum of the states: + cout << "Energy and momentum of states:" << endl; + cout << "gstate.E = " << gstate.E << "\testate.E = " << estate.E << endl; + cout << "Excitation energy = estate.E - gstate.E = " << estate.E - gstate.E << endl; + cout << "Excitation momentum (mod 2 pi) = estate.K - gstate.K = " << estate.K - gstate.K << endl; + + + // Computing matrix elements: CAREFUL: magnetizations must be consistent with operator (error is flagged). + cout << "Matrix elements: " << endl; + // The logarithm of the matrix elements are computed as complex numbers: + cout << "ln_Sz (gstate, estate) = " << ln_Sz_ME (gstate, estate) << endl; + + // The other possible matrix element call is for the Smin matrix element, + cout << "ln_Smin (gstate, estate2) = " << ln_Smin_ME (gstate, estate2) << endl; + + + clock_t StopTime = clock(); + + //cout << "Total time: " << (StopTime - StartTime)/*/CLOCKS_PER_SEC*/ << " hundreths of a second." + cout << "Total time: " << DP(StopTime - StartTime)/CLOCKS_PER_SEC << " seconds." + << endl; + + return(0); +} diff --git a/ABACUS_Usage_Example_LiebLin.cc b/ABACUS_Usage_Example_LiebLin.cc new file mode 100644 index 0000000..fc8447e --- /dev/null +++ b/ABACUS_Usage_Example_LiebLin.cc @@ -0,0 +1,108 @@ +/********************************************************** + +This software is part of J.-S. Caux's ABACUS library. + +Copyright (c) J.-S. Caux. + +----------------------------------------------------------- + +File: ABACUS_Usage_Example_LiebLin.cc + +Purpose: examples of calculations for Lieb-Liniger + +***********************************************************/ + +#include "ABACUS.h" + +using namespace std; +using namespace ABACUS; + + +int main() +{ + clock_t StartTime = clock(); + + DP c_int = 2.0; + DP L = 3.0; + int N = 3; + DP nbar_required = 1.0; + DP kBT = 4.0; + int Npts = 4*N; + DP req_diff = 1.0e-4; + int Max_Secs = 60; + + + if (true) { // State-by-state checks + + DP c_int = 4.0; + DP L = 16.0; + int N = 16; + + LiebLin_Bethe_State gstate (c_int, L, N); + gstate.Compute_All(true); + cout << setprecision(16) << gstate << endl; + + LiebLin_Bethe_State estate (c_int, L, N-1); + //estate.Set_to_ids (1LL, 1LL, 2LL, 0LL); + //estate.Set_to_Label ("64_2_028ysn1", gstate.Ix2); + //for (int i = 0; i < N; ++i) estate.Ix2[i] += 2; + // estate.Ix2[0] -= 8; + // estate.Ix2[1] -= 4; + // estate.Ix2[N-3] += 2; + // estate.Ix2[N-2] += 4; + // estate.Ix2[N-1] += 6; + // estate.Set_Label_from_Ix2(gstate.Ix2); + estate.Compute_All(true); + cout << setprecision(16) << estate << endl; + + //Scan_LiebLin ('o', estate, "28_3_i3y55yf3", -2*N, 2*N, 60, 1.0e+6, 0, 0, 1); + //stringstream filenameprefix; + //Data_File_Name (filenameprefix, 'd', -2*N, 2*N, 0.0, estate, estate, "28_3_i3y55yf3"); + //string prefix = filenameprefix.str(); + + DP ommin = 0.0; DP ommax = 10.0; DP gwidth = 1.0;// meaningless + int Nom = 10; + DP normalization = twoPI * L; + + //Smoothen_RAW_into_SF (filenameprefix.str(), -2*N, 2*N, 0, ommin, ommax, Nom, gwidth, normalization, L); + + //Density-density matrix elements: + // cout << setprecision(16) << "omega = " << estate.E - gstate.E << "\t" + // << exp(real(ln_Density_ME(gstate, estate))) << "\t" << exp(real(ln_Density_ME(estate, gstate))) << endl; + + //Field operator matrix elements: + cout << "omega\tiK\t< estate | Psi | gstate > matrix element:" << endl; + cout << setprecision(16) << estate.E - gstate.E << "\t" << estate.iK - gstate.iK << "\t" + << exp(ln_Psi_ME(estate, gstate)) << endl; + + //LiebLin_Bethe_State flipstate = estate; + //flipstate.Parity_Flip(); + //cout << "Flipping: " << endl; + //cout << "omega = " << flipstate.E - gstate.E << "\t" << exp(real(ln_Density_ME(gstate, flipstate))) << endl; + + // Finite T checks: + DP kBT = 1.0; + + // Delta is the number of sites involved in the smoothing of the entropy + //int Delta = int(sqrt(N))/2;//6;//N/20; + //DP epsilon = log(L)/L; + DP epsilon = log(L)/L; + + // Construct the finite-size saddle-point state: + //LiebLin_Bethe_State spstate = Canonical_Saddle_Point_State (c_int, L, N, kBT, Delta); + //LiebLin_Bethe_State spstate = Canonical_Saddle_Point_State (c_int, L, N, kBT, epsilon); + //spstate.Compute_All(true); + + //cout << spstate << endl; + + } + + + clock_t StopTime = clock(); + + //cout << "Total time: " << (StopTime - StartTime)/*/CLOCKS_PER_SEC*/ << " hundreths of a second." + cout << "Total time: " << DP(StopTime - StartTime)/CLOCKS_PER_SEC << " seconds." + << endl; + + return(0); +} diff --git a/include/ABACUS_LiebLin.h b/include/ABACUS_LiebLin.h index 349f8bc..03cda15 100644 --- a/include/ABACUS_LiebLin.h +++ b/include/ABACUS_LiebLin.h @@ -29,6 +29,14 @@ namespace ABACUS { //*********************************************************************** + /** + Eigenstates of the Lieb-Liniger model in the repulsive \f$(c > 0)\f$ regime. + + For numerical convenience, rapidities are rescaled by the interaction parameter \f$c\f$. + In equations throughout this documentation, \f$\tilde{\lambda} \equiv \lambda/c\f$. + In the code, these rescaled rapidities are denoted `lambdaoc` (read: "lambda over c"). + */ + class LiebLin_Bethe_State { public: @@ -104,7 +112,8 @@ namespace ABACUS { }; inline bool Is_Inner_Skeleton (LiebLin_Bethe_State& State) { - return (State.N >= 2 && (State.Ix2[0] == State.Ix2[1] - 2 || State.Ix2[State.N-1] == State.Ix2[State.N-2] + 2)); + return (State.N >= 2 && (State.Ix2[0] == State.Ix2[1] - 2 + || State.Ix2[State.N-1] == State.Ix2[State.N-2] + 2)); }; inline bool Is_Outer_Skeleton (LiebLin_Bethe_State& State) { return (State.N >= 2 && State.Ix2[0] == LIEBLIN_Ix2_MIN + (State.N % 2) + 1 @@ -115,12 +124,14 @@ namespace ABACUS { && State.Ix2[State.N-1] == LIEBLIN_Ix2_MAX - (State.N % 2) - 1); }; - inline bool Force_Descent (char whichDSF, LiebLin_Bethe_State& ScanState, LiebLin_Bethe_State& RefState, + inline bool Force_Descent (char whichDSF, LiebLin_Bethe_State& ScanState, + LiebLin_Bethe_State& RefState, int desc_type_required, int iKmod, DP Chem_Pot) { bool forcedesc = false; - // Force descent if we're computing density-density, we're at zero momentum and we're descending with momentum preserved: + // Force descent if we're computing density-density, we're at zero momentum + // and we're descending with momentum preserved: if (whichDSF == 'd' && ScanState.iK == RefState.iK && desc_type_required > 8) forcedesc = true; // For BEC to c > 0 quench, g2(x=0): force first step @@ -135,7 +146,8 @@ namespace ABACUS { // FUNCTION DECLARATIONS: DP Chemical_Potential (LiebLin_Bethe_State& RefState); - DP Sumrule_Factor (char whichDSF, LiebLin_Bethe_State& RefState, DP Chem_Pot, int iKmin, int iKmax); + DP Sumrule_Factor (char whichDSF, LiebLin_Bethe_State& RefState, DP Chem_Pot, + int iKmin, int iKmax); void Evaluate_F_Sumrule (char whichDSF, const LiebLin_Bethe_State& RefState, DP Chem_Pot, int iKmin, int iKmax, const char* RAW_Cstr, const char* FSR_Cstr); void Evaluate_F_Sumrule (std::string prefix, char whichDSF, const LiebLin_Bethe_State& RefState, @@ -157,14 +169,23 @@ namespace ABACUS { std::complex ln_Psi_ME (LiebLin_Bethe_State& lstate, LiebLin_Bethe_State& rstate); std::complex ln_g2_ME (LiebLin_Bethe_State& mu, LiebLin_Bethe_State& lambda); - DP Compute_Matrix_Element_Contrib (char whichDSF, int iKmin, int iKmax, LiebLin_Bethe_State& LeftState, - LiebLin_Bethe_State& RefState, DP Chem_Pot, std::stringstream& DAT_outfile); + DP Compute_Matrix_Element_Contrib (char whichDSF, int iKmin, int iKmax, + LiebLin_Bethe_State& LeftState, + LiebLin_Bethe_State& RefState, DP Chem_Pot, + std::stringstream& DAT_outfile); DP LiebLin_Twisted_lnnorm (Vect >& lambdaoc, double cxL); - std::complex LiebLin_Twisted_ln_Overlap (DP expbeta, Vect lstate_lambdaoc, DP lstate_lnnorm, LiebLin_Bethe_State& rstate); - std::complex LiebLin_Twisted_ln_Overlap (std::complex expbeta, Vect > lstate_lambdaoc, DP lstate_lnnorm, LiebLin_Bethe_State& rstate); - std::complex LiebLin_ln_Overlap (Vect lstate_lambdaoc, DP lstate_lnnorm, LiebLin_Bethe_State& rstate); - std::complex LiebLin_ln_Overlap (Vect > lstate_lambdaoc, DP lstate_lnnorm, LiebLin_Bethe_State& rstate); + std::complex LiebLin_Twisted_ln_Overlap (DP expbeta, Vect lstate_lambdaoc, + DP lstate_lnnorm, + LiebLin_Bethe_State& rstate); + std::complex LiebLin_Twisted_ln_Overlap (std::complex expbeta, + Vect > lstate_lambdaoc, + DP lstate_lnnorm, + LiebLin_Bethe_State& rstate); + std::complex LiebLin_ln_Overlap (Vect lstate_lambdaoc, DP lstate_lnnorm, + LiebLin_Bethe_State& rstate); + std::complex LiebLin_ln_Overlap (Vect > lstate_lambdaoc, DP lstate_lnnorm, + LiebLin_Bethe_State& rstate); // In src/LiebLin_Tgt0.cc: DP Entropy (LiebLin_Bethe_State& RefState); diff --git a/src/LIEBLIN/LiebLin_Bethe_State.cc b/src/LIEBLIN/LiebLin_Bethe_State.cc index 7805f59..088f559 100644 --- a/src/LIEBLIN/LiebLin_Bethe_State.cc +++ b/src/LIEBLIN/LiebLin_Bethe_State.cc @@ -24,18 +24,18 @@ namespace ABACUS { LiebLin_Bethe_State::LiebLin_Bethe_State () : c_int (0.0), L(0.0), cxL(0.0), N(0), - Ix2_available(Vect(0, 1)), index_first_hole_to_right (Vect(0,1)), displacement (Vect(0,1)), - Ix2(Vect(0, 1)), lambdaoc(Vect(0.0, 1)), - S(Vect(0.0, 1)), dSdlambdaoc(Vect(0.0, 1)), - diffsq(0.0), prec(ITER_REQ_PREC_LIEBLIN), conv(0), iter_Newton(0), E(0.0), iK(0), K(0.0), lnnorm(-100.0) + Ix2_available(Vect(0, 1)), index_first_hole_to_right (Vect(0,1)), + displacement (Vect(0,1)), Ix2(Vect(0, 1)), lambdaoc(Vect(0.0, 1)), + S(Vect(0.0, 1)), dSdlambdaoc(Vect(0.0, 1)), diffsq(0.0), prec(ITER_REQ_PREC_LIEBLIN), + conv(0), iter_Newton(0), E(0.0), iK(0), K(0.0), lnnorm(-100.0) { stringstream Nout; Nout << N; label = Nout.str() + LABELSEP + ABACUScoding[0] + LABELSEP; } LiebLin_Bethe_State::LiebLin_Bethe_State (DP c_int_ref, DP L_ref, int N_ref) : c_int(c_int_ref), L(L_ref), cxL(c_int_ref * L_ref), N(N_ref), - Ix2_available(Vect(0, 2)), index_first_hole_to_right (Vect(0,N)), displacement (Vect(0,N)), - Ix2(Vect(0, N)), lambdaoc(Vect(0.0, N)), + Ix2_available(Vect(0, 2)), index_first_hole_to_right (Vect(0,N)), + displacement (Vect(0,N)), Ix2(Vect(0, N)), lambdaoc(Vect(0.0, N)), S(Vect(0.0, N)), dSdlambdaoc(Vect(0.0, N)), diffsq(0.0), prec(ABACUS::max(1.0, 1.0/(c_int * c_int)) * ITER_REQ_PREC_LIEBLIN), conv(0), iter_Newton(0), E(0.0), iK(0), K(0.0), lnnorm(-100.0) @@ -93,7 +93,8 @@ namespace ABACUS { if (N != OriginStateIx2.size()) { cout << label_ref << endl; cout << labeldata.M << endl; - ABACUSerror("Trying to set an incorrect label on LiebLin_Bethe_State: N != OriginStateIx2.size()."); + ABACUSerror("Trying to set an incorrect label on LiebLin_Bethe_State: " + "N != OriginStateIx2.size()."); } label = label_ref; @@ -106,7 +107,8 @@ namespace ABACUS { // Now set the excitations: for (int iexc = 0; iexc < labeldata.nexc[0]; ++iexc) - for (int i = 0; i < N; ++i) if (Ix2[i] == labeldata.Ix2old[0][iexc]) Ix2[i] = labeldata.Ix2exc[0][iexc]; + for (int i = 0; i < N; ++i) + if (Ix2[i] == labeldata.Ix2old[0][iexc]) Ix2[i] = labeldata.Ix2exc[0][iexc]; // Now reorder the Ix2 to follow convention: Ix2.QuickSort(); @@ -143,10 +145,13 @@ namespace ABACUS { Ix2old_ref[0] = Vect(ABACUS::max(nexc_ref[0],1)); Ix2exc_ref[0] = Vect(ABACUS::max(nexc_ref[0],1)); int nexccheck = 0; - for (int i = 0; i < N; ++i) if (!OriginStateIx2.includes(Ix2[i])) Ix2exc_ref[0][nexccheck++] = Ix2[i]; - if (nexccheck != nexc_ref[0]) ABACUSerror("Counting excitations wrong (1) in LiebLin_Bethe_State::Set_Label_from_Ix2"); + for (int i = 0; i < N; ++i) + if (!OriginStateIx2.includes(Ix2[i])) Ix2exc_ref[0][nexccheck++] = Ix2[i]; + if (nexccheck != nexc_ref[0]) + ABACUSerror("Counting excitations wrong (1) in LiebLin_Bethe_State::Set_Label_from_Ix2"); nexccheck = 0; - for (int i = 0; i < N; ++i) if (!Ix2.includes (OriginStateIx2[i])) Ix2old_ref[0][nexccheck++] = OriginStateIx2[i]; + for (int i = 0; i < N; ++i) + if (!Ix2.includes (OriginStateIx2[i])) Ix2old_ref[0][nexccheck++] = OriginStateIx2[i]; if (nexccheck != nexc_ref[0]) { cout << "nexc_ref[0] = " << nexc_ref[0] << "\tnexccheck = " << nexccheck << endl; cout << OriginStateIx2 << endl; @@ -167,7 +172,8 @@ namespace ABACUS { void LiebLin_Bethe_State::Set_Label_Internals_from_Ix2 (const Vect& OriginStateIx2) { - if (N != OriginStateIx2.size()) ABACUSerror("N != OriginStateIx2.size() in Set_Label_Internals_from_Ix2."); + if (N != OriginStateIx2.size()) + ABACUSerror("N != OriginStateIx2.size() in Set_Label_Internals_from_Ix2."); Vect OriginStateIx2ordered = OriginStateIx2; OriginStateIx2ordered.QuickSort(); @@ -207,7 +213,8 @@ namespace ABACUS { for (int j = 0; j < N; ++j) { int i = 0; while (Ix2_available[i] < OriginStateIx2ordered[j]) i++; - // We now have Ix2_available[i] == OriginStateIx2[j]. Shift all Ix2_available to the right of this by 2; + // We now have Ix2_available[i] == OriginStateIx2[j]. + // Shift all Ix2_available to the right of this by 2; for (int i1 = i; i1 < navailable; ++i1) Ix2_available[i1] += 2; index_first_hole_to_right[j] = i; } @@ -218,10 +225,12 @@ namespace ABACUS { // Set displacement vector from the Ix2: for (int j = 0; j < N; ++j) { if (Ix2[j] < OriginStateIx2ordered[j]) { - // Ix2[j] must be equal to some OriginState_Ix2_available[i] for i < OriginState_index_first_hole_to_right[j] + // Ix2[j] must be equal to some OriginState_Ix2_available[i] + // for i < OriginState_index_first_hole_to_right[j] while (Ix2[j] != Ix2_available[index_first_hole_to_right[j] + displacement[j] ]) { if (index_first_hole_to_right[j] + displacement[j] == 0) { - cout << label << endl << j << endl << OriginStateIx2 << endl << Ix2 << endl << Ix2_available + cout << label << endl << j << endl << OriginStateIx2 << endl + << Ix2 << endl << Ix2_available << endl << index_first_hole_to_right << endl << displacement << endl; ABACUSerror("Going down too far in Set_Label_Internals..."); } @@ -233,7 +242,8 @@ namespace ABACUS { displacement[j] = 1; // start with this value to prevent segfault while (Ix2[j] != Ix2_available[index_first_hole_to_right[j] - 1 + displacement[j] ]) { if (index_first_hole_to_right[j] + displacement[j] == Ix2_available.size() - 1) { - cout << label << endl << j << endl << OriginStateIx2 << endl << Ix2 << endl << Ix2_available + cout << label << endl << j << endl << OriginStateIx2 << endl + << Ix2 << endl << Ix2_available << endl << index_first_hole_to_right << endl << displacement << endl; ABACUSerror("Going up too far in Set_Label_Internals..."); } @@ -245,7 +255,8 @@ namespace ABACUS { bool LiebLin_Bethe_State::Check_Admissibility (char whichDSF) { - //if (Ix2.min() < -13 || Ix2.max() > 13) return(false); // For testing with restricted Hilbert space + // For testing with restricted Hilbert space: + //if (Ix2.min() < -13 || Ix2.max() > 13) return(false); return(true); } @@ -253,7 +264,8 @@ namespace ABACUS { { // This function finds the rapidities of the eigenstate - lnnorm = -100.0; // sentinel value, recalculated if Newton method used in the last step of iteration. + // sentinel value, recalculated if Newton method used in the last step of iteration. + lnnorm = -100.0; diffsq = 1.0; @@ -289,15 +301,20 @@ namespace ABACUS { return nonan; } + /** + For the repulsive Lieb-Liniger model, all eigenstates have real rapidities + (there are no strings). All string deviations are thus set to zero. + */ DP LiebLin_Bethe_State::String_delta() { - return(0.0); // no strings (thus no deviations) in replusive LiebLin + return(0.0); } + /** + Checks whether the quantum numbers are symmetrically distributed: \f$\{ I \} = \{ -I \}\f$. + */ bool LiebLin_Bethe_State::Check_Symmetry () { - // Checks whether the I's are symmetrically distributed. - bool symmetric_state = true; Vect Ix2check = Ix2; @@ -309,6 +326,18 @@ namespace ABACUS { return(symmetric_state); } + /** + This function computes the log of the norm of a `LiebLin_Bethe_State` instance. + This is obtained from the Gaudin-Korepin norm formula + \f{eqnarray*}{ + \langle 0 | \prod_{j=0}^{N-1} C(\lambda_j) \prod_{j=0}^{N-1} B(\lambda_j) | 0 \rangle + = \prod_{j=0}^{N-2}\prod_{k=j+1}^{N-1} + \frac{(\lambda_j - \lambda_k )^2 + c^2}{(\lambda_j - \lambda_k)^2} + \times ~\mbox{Det}_N G + \f} + in which \f$G\f$ is the Gaudin matrix computed by + LiebLin_Bethe_State::Build_Reduced_Gaudin_Matrix. + */ void LiebLin_Bethe_State::Compute_lnnorm () { if (lnnorm == -100.0) { // else Gaudin part already calculated by Newton method @@ -317,7 +346,13 @@ namespace ABACUS { (*this).Build_Reduced_Gaudin_Matrix(Gaudin_Red); - lnnorm = real(lndet_LU_dstry(Gaudin_Red)); + complex lnnorm_CX = lndet_LU_dstry(Gaudin_Red); + if (abs(imag(lnnorm_CX)) > 1.0e-6) { + cout << "lnnorm_CX = " << lnnorm_CX << endl; + ABACUSerror("Gaudin norm of LiebLin_Bethe_State should be real and positive."); + } + + lnnorm = real(lnnorm_CX); // Add the pieces outside of Gaudin determinant @@ -329,7 +364,15 @@ namespace ABACUS { return; } - void LiebLin_Bethe_State::Compute_All (bool reset_rapidities) // solves BAE, computes E, K and lnnorm + /** + This function solves the Bethe equations, computes the energy, momentum and state norm. + It calls LiebLin_Bethe_State::Find_Rapidities, LiebLin_Bethe_State::Compute_Energy, + LiebLin_Bethe_State::Compute_Momentum and LiebLin_Bethe_State::Compute_lnnorm. + + Assumptions: + - the set of quantum numbers is admissible. + */ + void LiebLin_Bethe_State::Compute_All (bool reset_rapidities) { (*this).Find_Rapidities (reset_rapidities); if (conv == 1) { @@ -345,7 +388,8 @@ namespace ABACUS { if (cxL >= 1.0) for (int a = 0; a < N; ++a) lambdaoc[a] = PI * Ix2[a]/cxL; - // For small values of c, use better approximation using approximate zeroes of Hermite polynomials: see Gaudin eqn 4.71. + // For small values of c, use better approximation using approximate + // zeroes of Hermite polynomials: see Gaudin eqn 4.71. if (cxL < 1.0) { DP oneoversqrtcLN = 1.0/pow(cxL * N, 0.5); for (int a = 0; a < N; ++a) lambdaoc[a] = oneoversqrtcLN * PI * Ix2[a]; @@ -403,7 +447,8 @@ namespace ABACUS { dSdlambdaoc[j] += 1.0/((lambdaoc[j] - lambdaoc[k]) * (lambdaoc[j] - lambdaoc[k]) + 1.0); dSdlambdaoc[j] *= 2.0/(PI * cxL); - dlambdaoc[j] = (PI*Ix2[j]/cxL - S[j] + lambdaoc[j] * dSdlambdaoc[j])/(1.0 + dSdlambdaoc[j]) - lambdaoc[j]; + dlambdaoc[j] = (PI*Ix2[j]/cxL - S[j] + lambdaoc[j] * dSdlambdaoc[j]) + /(1.0 + dSdlambdaoc[j]) - lambdaoc[j]; } @@ -421,16 +466,18 @@ namespace ABACUS { return; } + /** + Performs one step of the matrix Newton method on the rapidities. + */ void LiebLin_Bethe_State::Iterate_BAE_Newton (DP damping) { - // does one step of a Newton method on the rapidities... - Vect_DP RHSBAE (0.0, N); // contains RHS of BAEs Vect_DP dlambdaoc (0.0, N); // contains delta lambdaoc computed from Newton's method SQMat_DP Gaudin (0.0, N); Vect_INT indx (N); DP sumtheta = 0.0; - int atanintshift = 0; // for large |lambda|, use atan (lambda) = sgn(lambda) pi/2 - atan(1/lambda) + // for large |lambda|, use atan (lambda) = sgn(lambda) pi/2 - atan(1/lambda) + int atanintshift = 0; DP lambdahere = 0.0; // Compute the RHS of the BAEs: @@ -463,11 +510,13 @@ namespace ABACUS { lubksb (Gaudin, indx, dlambdaoc); bool ordering_changed = false; - for (int j = 0; j < N-1; ++j) if (lambdaoc[j] + dlambdaoc[j] > lambdaoc[j+1] + dlambdaoc[j+1]) ordering_changed = true; + for (int j = 0; j < N-1; ++j) + if (lambdaoc[j] + dlambdaoc[j] > lambdaoc[j+1] + dlambdaoc[j+1]) ordering_changed = true; // To prevent Newton from diverging, we limit the size of the rapidity changes. // The leftmost and rightmost rapidities can grow by one order of magnitude per iteration step. - if (ordering_changed) { // We explicitly ensure that the ordering remains correct after the iteration step. + if (ordering_changed) { + // We explicitly ensure that the ordering remains correct after the iteration step. bool ordering_still_changed = false; DP maxdlambdaoc = 0.0; do { @@ -552,16 +601,38 @@ namespace ABACUS { K = 2.0 * iK * PI/L; } + /** + The fundamental scattering kernel for Lieb-Liniger, + \f[ + K (\tilde{\lambda}) = \frac{ 2 }{\tilde{\lambda}^2 + 1} + \f] + given two indices for the rapidities as arguments. + */ DP LiebLin_Bethe_State::Kernel (int a, int b) { return(2.0/(pow(lambdaoc[a] - lambdaoc[b], 2.0) + 1.0)); } + /** + The fundamental scattering kernel for Lieb-Liniger, + \f[ + K (\tilde{\lambda}) = \frac{ 2 }{\tilde{\lambda}^2 + 1} + \f] + given the rapidity difference as argument. + */ DP LiebLin_Bethe_State::Kernel (DP lambdaoc_ref) { return(2.0/(lambdaoc_ref * lambdaoc_ref + 1.0)); } + /** + This function constructs the Gaudin matrix \f$G\f$ + which is defined as the Hessian of the Yang-Yang action \f$S^{YY}\f$, + \f[ + G_{jk} = \delta_{jk} \left\{ cL + \sum_{l} K (\tilde{\lambda}_j - \tilde{\lambda}_l) \right\} + - K (\tilde{\lambda}_j - \tilde{\lambda}_k). + \f] + */ void LiebLin_Bethe_State::Build_Reduced_Gaudin_Matrix (SQMat& Gaudin_Red) { @@ -575,7 +646,8 @@ namespace ABACUS { if (j == k) { sum_Kernel = 0.0; - for (int kp = 0; kp < N; ++kp) if (j != kp) sum_Kernel += Kernel (lambdaoc[j] - lambdaoc[kp]); + for (int kp = 0; kp < N; ++kp) + if (j != kp) sum_Kernel += Kernel (lambdaoc[j] - lambdaoc[kp]); Gaudin_Red[j][k] = cxL + sum_Kernel; } @@ -586,6 +658,11 @@ namespace ABACUS { return; } + /** + For a summetric `LiebLin_Bethe_State`, this function computes the Gaudin matrix + as an \f$N/2 \times N/2\f$ matrix to accelerate computations. + */ + void LiebLin_Bethe_State::Build_Reduced_BEC_Quench_Gaudin_Matrix (SQMat& Gaudin_Red) { // Passing a matrix of dimension N/2 @@ -612,7 +689,8 @@ namespace ABACUS { sum_Kernel = 0.0; for (int kp = N/2; kp < N; ++kp) if (j + N/2 != kp) - sum_Kernel += Kernel (lambdaoc[j+N/2] - lambdaoc[kp]) + Kernel (lambdaoc[j+N/2] + lambdaoc[kp]); + sum_Kernel += Kernel (lambdaoc[j+N/2] - lambdaoc[kp]) + + Kernel (lambdaoc[j+N/2] + lambdaoc[kp]); Gaudin_Red[j][k] = cxL + sum_Kernel; } @@ -625,10 +703,13 @@ namespace ABACUS { } - void LiebLin_Bethe_State::Annihilate_ph_pair (int ipart, int ihole, const Vect& OriginStateIx2) + /** + This function changes the Ix2 of a given state by annihilating a particle and hole + pair specified by ipart and ihole (counting from the left, starting with index 0). + */ + void LiebLin_Bethe_State::Annihilate_ph_pair (int ipart, int ihole, + const Vect& OriginStateIx2) { - // This function changes the Ix2 of a given state by annihilating a particle and hole - // pair specified by ipart and ihole (counting from the left, starting with index 0). State_Label_Data currentdata = Read_State_Label ((*this).label, OriginStateIx2); @@ -645,8 +726,10 @@ namespace ABACUS { int ntypespresent = 1; // only one type for LiebLin Vect > Ix2old_new(ntypespresent); Vect > Ix2exc_new(ntypespresent); - for (int it = 0; it < ntypespresent; ++it) Ix2old_new[it] = Vect(ABACUS::max(nexc_new[it],1)); - for (int it = 0; it < ntypespresent; ++it) Ix2exc_new[it] = Vect(ABACUS::max(nexc_new[it],1)); + for (int it = 0; it < ntypespresent; ++it) + Ix2old_new[it] = Vect(ABACUS::max(nexc_new[it],1)); + for (int it = 0; it < ntypespresent; ++it) + Ix2exc_new[it] = Vect(ABACUS::max(nexc_new[it],1)); // Copy earlier data in, leaving out ipart and ihole: for (int it = 0; it < ntypespresent; ++it) { @@ -661,9 +744,12 @@ namespace ABACUS { (*this).Set_to_Label (Return_State_Label(newdata, OriginStateIx2)); } + /** + Performs a spatial inversion of a `LiebLin_Bethe_State`, mapping all quantum numbers + and rapidities to minus themselves. + */ void LiebLin_Bethe_State::Parity_Flip () { - // For simplicity, we don't redo base_id, type_id, id. Vect_INT Ix2buff = Ix2; Vect_DP lambdaocbuff = lambdaoc; for (int i = 0; i < N; ++i) Ix2[i] = -Ix2buff[N - 1 - i]; @@ -672,16 +758,20 @@ namespace ABACUS { K = -K; } + /** + Send information about a `LiebLin_Bethe_State` to the output stream. + */ std::ostream& operator<< (std::ostream& s, const LiebLin_Bethe_State& state) { - s << endl << "******** State for c = " << state.c_int << " L = " << state.L << " N = " << state.N - << " with label " << state.label << " ********" << endl; + s << endl << "******** State for c = " << state.c_int << " L = " << state.L + << " N = " << state.N << " with label " << state.label << " ********" << endl; s << "Ix2:" << endl; for (int j = 0; j < state.N; ++j) s << state.Ix2[j] << " "; s << endl << "lambdaocs:" << endl; for (int j = 0; j < state.N; ++j) s << state.lambdaoc[j] << " "; s << endl << "conv = " << state.conv << " iter_Newton = " << state.iter_Newton << endl; - s << "E = " << state.E << " iK = " << state.iK << " K = " << state.K << " lnnorm = " << state.lnnorm << endl; + s << "E = " << state.E << " iK = " << state.iK << " K = " << state.K + << " lnnorm = " << state.lnnorm << endl; return(s); } diff --git a/src/LIEBLIN/ln_Density_ME.cc b/src/LIEBLIN/ln_Density_ME.cc index 3e91d81..32ff25e 100644 --- a/src/LIEBLIN/ln_Density_ME.cc +++ b/src/LIEBLIN/ln_Density_ME.cc @@ -19,35 +19,51 @@ using namespace ABACUS; namespace ABACUS { + /** + Given a index \f$j\f$ together with a bra and ket respectively parametrized + by sets of rapidities \f$\{ \mu \}\f$ and \f$\{ \lambda \}\f$, + this function returns the product (see e.g. + [J.-S. Caux et al, JSTAT P01008 (2007)](https://doi.org/10.1088/1742-5468/2007/01/P01008)) + \f[ + V_j^+ \equiv \prod_l + \frac{\tilde{\mu}_l - \tilde{\lambda}_j + i}{\tilde{\lambda}_l - \tilde{\lambda}_j + i}. + \f] + */ complex Fn_V (int j, LiebLin_Bethe_State& lstate, LiebLin_Bethe_State& rstate) { complex result = 1.0; for (int m = 0; m < lstate.N; ++m) { - result *= (lstate.lambdaoc[m] - rstate.lambdaoc[j] + II)/(rstate.lambdaoc[m] - rstate.lambdaoc[j] + II); + result *= (lstate.lambdaoc[m] - rstate.lambdaoc[j] + II) + /(rstate.lambdaoc[m] - rstate.lambdaoc[j] + II); } return(result); } + /** + Computes the log of the density operator \f$\hat{\rho}(x = 0)\f$ matrix element + between lstate (bra) and rstate (ket). If lstate == rstate, density matrix element = N/L; + if momentum difference is zero but states are different, then matrix element is zero + (but this function then returns the log artificially set to -200.0). + */ complex ln_Density_ME (LiebLin_Bethe_State& lstate, LiebLin_Bethe_State& rstate) { - // Computes the log of the density operator \rho(x = 0) matrix element between lstate and rstate. - - // If we have lstate == rstate, density matrix element = N/L: if (lstate.Ix2 == rstate.Ix2) return(log(lstate.N/lstate.L)); - // If momentum difference is zero but states are different, then matrix element is zero: else if (lstate.iK == rstate.iK) return(-200.0); // so exp(.) is zero SQMat_DP one_plus_U (0.0, lstate.N); Vect_CX Vplus (0.0, lstate.N); // contains V^+_j; V^-_j is the conjugate - Vect_CX ln_Fn_Prod (0.0, lstate.N); // product_{m\neq j} (\mu_m - \lambdaoc_j)/(\lambdaoc_m - \lambdaoc_j) + // prod_{m\neq j} (\mu_m - \lambdaoc_j)/(\lambdaoc_m - \lambdaoc_j): + Vect_CX ln_Fn_Prod (0.0, lstate.N); Vect_DP rKern (0.0, lstate.N); // K(lambdaoc_j - lambdaoc_(p == arbitrary)) + // "Phantom" rapidity in ME expression. Choice doesn't matter, see 1990_Slavnov_TMP_82 + // after (3.8). Choose rapidity around the middle. //int p = 0; - int p = rstate.N/2-1; // choice doesn't matter, see 1990_Slavnov_TMP_82 after (3.8). Choose rapidity around the middle. + int p = rstate.N/2-1; DP Kout = lstate.K - rstate.K; @@ -62,8 +78,9 @@ namespace ABACUS { for (int a = 0; a < lstate.N; ++a) for (int b = 0; b < lstate.N; ++b) - one_plus_U[a][b] = (a == b ? 1.0 : 0.0) + 0.5 * ((lstate.lambdaoc[a] - rstate.lambdaoc[a])/imag(Vplus[a])) - * real(exp(ln_Fn_Prod[a])) * (rstate.Kernel(a,b) - rKern[b]); // BUGRISK: why real here? + one_plus_U[a][b] = (a == b ? 1.0 : 0.0) + + 0.5 * ((lstate.lambdaoc[a] - rstate.lambdaoc[a])/imag(Vplus[a])) + * real(exp(ln_Fn_Prod[a])) * (rstate.Kernel(a,b) - rKern[b]); complex ln_ddalpha_sigma = lndet_LU_dstry(one_plus_U); diff --git a/src/LIEBLIN/ln_Psi_ME.cc b/src/LIEBLIN/ln_Psi_ME.cc index 36d4ae3..38c11cb 100644 --- a/src/LIEBLIN/ln_Psi_ME.cc +++ b/src/LIEBLIN/ln_Psi_ME.cc @@ -19,62 +19,78 @@ using namespace ABACUS; namespace ABACUS { + /** + Given a index \f$j\f$ together with a bra and ket respectively parametrized + by sets of rapidities \f$\{ \mu \}\f$ and \f$\{ \lambda \}\f$, + this function returns the product (see e.g. + [J.-S. Caux et al, JSTAT P01008 (2007)](https://doi.org/10.1088/1742-5468/2007/01/P01008)) + \f[ + V_j^+ \equiv \prod_l + \frac{\tilde{\mu}_l - \tilde{\lambda}_j + i}{\tilde{\lambda}_l - \tilde{\lambda}_j + i}. + \f] + */ complex Fn_V_Psi (int j, LiebLin_Bethe_State& lstate, LiebLin_Bethe_State& rstate) { // lstate has one particle less than rstate - complex result_num = 1.0; - complex result_den = 1.0; + complex result = 1.0; for (int m = 0; m < lstate.N; ++m) - result_num *= (lstate.lambdaoc[m] - rstate.lambdaoc[j] + II); + result *= (lstate.lambdaoc[m] - rstate.lambdaoc[j] + II) + /(rstate.lambdaoc[m] - rstate.lambdaoc[j] + II); + result /= (rstate.lambdaoc[rstate.N - 1] - rstate.lambdaoc[j] + II); - for (int m = 0; m < rstate.N; ++m) - result_den *= (rstate.lambdaoc[m] - rstate.lambdaoc[j] + II); - - return(result_num/result_den); + return(result); } + /** + Computes the log of the annihilation operator matrix element between lstate and rstate. + */ complex ln_Psi_ME (LiebLin_Bethe_State& lstate, LiebLin_Bethe_State& rstate) { - // Computes the log of the annihilation operator matrix element between lstate and rstate. - - if (lstate.N + 1 != rstate.N) ABACUSerror("Wrong particle numbers in left and right states for Psi FF."); + if (lstate.N + 1 != rstate.N) + ABACUSerror("Wrong particle numbers in left and right states for Psi FF."); SQMat_DP U_Psi (0.0, lstate.N); Vect_CX Vplus (0.0, lstate.N); // contains V^+_j; V^-_j is the conjugate - Vect_DP Fn_Prod (0.0, lstate.N); // product_{m} (\mu_m - \lambdaoc_j)/product_{m neq j} (\lambdaoc_m - \lambdaoc_j) + // product_{m} (\mu_m - \lambdaoc_j)/product_{m neq j} (\lambdaoc_m - \lambdaoc_j) + Vect_DP Fn_Prod (0.0, lstate.N); Vect_DP rKern (0.0, lstate.N); // K(lambdaoc_j - lambdaoc_(p == N)) int p = rstate.N - 1; for (int a = 0; a < lstate.N; ++a) { Vplus[a] = Fn_V_Psi (a, lstate, rstate); - Fn_Prod[a] = (lstate.lambdaoc[a] - rstate.lambdaoc[a])/(rstate.lambdaoc[rstate.N - 1] - rstate.lambdaoc[a]); + Fn_Prod[a] = (lstate.lambdaoc[a] - rstate.lambdaoc[a]) + /(rstate.lambdaoc[rstate.N - 1] - rstate.lambdaoc[a]); for (int m = 0; m < lstate.N; ++m) - if (m != a) Fn_Prod[a] *= (lstate.lambdaoc[m] - rstate.lambdaoc[a])/(rstate.lambdaoc[m] - rstate.lambdaoc[a]); + if (m != a) Fn_Prod[a] *= (lstate.lambdaoc[m] - rstate.lambdaoc[a]) + /(rstate.lambdaoc[m] - rstate.lambdaoc[a]); rKern[a] = rstate.Kernel (a, p); } for (int a = 0; a < lstate.N; ++a) for (int b = 0; b < lstate.N; ++b) - U_Psi[a][b] = (a == b ? 2.0 * imag(Vplus[a]) : 0.0) + Fn_Prod[a] * (rstate.Kernel(a,b) - rKern[b]); + U_Psi[a][b] = (a == b ? 2.0 * imag(Vplus[a]) : 0.0) + + Fn_Prod[a] * (rstate.Kernel(a,b) - rKern[b]); complex ln_det_U_Psi = lndet_LU_dstry(U_Psi); complex ln_prod_lambdaocsq_plus_one = 0.0; for (int a = 0; a < rstate.N - 1; ++a) for (int b = a; b < rstate.N; ++b) - ln_prod_lambdaocsq_plus_one += log(complex((rstate.lambdaoc[a] - rstate.lambdaoc[b]) - * (rstate.lambdaoc[a] - rstate.lambdaoc[b])) + 1.0); + ln_prod_lambdaocsq_plus_one += + log(complex((rstate.lambdaoc[a] - rstate.lambdaoc[b]) + * (rstate.lambdaoc[a] - rstate.lambdaoc[b])) + 1.0); complex ln_prod_lambdaoca_min_mub = 0.0; for (int a = 0; a < rstate.N; ++a) for (int b = 0; b < lstate.N; ++b) ln_prod_lambdaoca_min_mub += log(complex(rstate.lambdaoc[a] - lstate.lambdaoc[b])); - return (ln_det_U_Psi + 0.5 * log(lstate.c_int) + ln_prod_lambdaocsq_plus_one - ln_prod_lambdaoca_min_mub + return (ln_det_U_Psi + 0.5 * log(lstate.c_int) + + ln_prod_lambdaocsq_plus_one - ln_prod_lambdaoca_min_mub - 0.5 * (lstate.lnnorm + rstate.lnnorm)); }