diff --git a/.clang-format b/.clang-format new file mode 100644 index 0000000..10b3840 --- /dev/null +++ b/.clang-format @@ -0,0 +1,45 @@ +BasedOnStyle: LLVM + +AccessModifierOffset: 0 +AlignAfterOpenBracket: Align +AlignConsecutiveAssignments: true +AlignConsecutiveDeclarations: false +AlignEscapedNewlinesLeft: false +AlignOperands: false +AlignTrailingComments: true +AllowAllParametersOfDeclarationOnNextLine: false +AllowShortBlocksOnASingleLine: true +AllowShortCaseLabelsOnASingleLine: true +AllowShortFunctionsOnASingleLine: All +AllowShortIfStatementsOnASingleLine: true +AllowShortLoopsOnASingleLine: true +AlwaysBreakBeforeMultilineStrings: true +AlwaysBreakTemplateDeclarations: false +BinPackArguments: true +BinPackParameters: true +BreakBeforeBinaryOperators: NonAssignment +BreakBeforeBraces: Attach +BreakBeforeTernaryOperators: false +BreakConstructorInitializersBeforeComma: false +BreakStringLiterals: false +ColumnLimit: 150 +ConstructorInitializerAllOnOneLineOrOnePerLine: true +ConstructorInitializerIndentWidth: 3 +ContinuationIndentWidth: 3 +Cpp11BracedListStyle: true +DerivePointerBinding : false +IndentCaseLabels: true +IndentWidth: 2 +Language: Cpp +MaxEmptyLinesToKeep: 1 +NamespaceIndentation : All +PointerAlignment: Right +ReflowComments: false +SortIncludes: false +SpaceAfterControlStatementKeyword: true +SpaceBeforeAssignmentOperators: true +SpaceInEmptyParentheses: false +SpacesInParentheses: false +Standard: Cpp11 +TabWidth: 2 +UseTab: Never diff --git a/CMakeLists.txt b/CMakeLists.txt new file mode 100644 index 0000000..5d7bd13 --- /dev/null +++ b/CMakeLists.txt @@ -0,0 +1,60 @@ +# Start configuration +cmake_minimum_required(VERSION 2.8.7) + +# Version number of the application +project(pomerol2triqs CXX) +set(POMEROL2TRIQS_VERSION "0.1") + +# Build the optimized version by default +if(NOT CMAKE_BUILD_TYPE) + set(CMAKE_BUILD_TYPE Release) +endif() + +# We use shared libraries +option(BUILD_SHARED_LIBS "Build shared libraries" ON) + +# Append TRIQS installed files to the cmake load path +list(APPEND CMAKE_MODULE_PATH ${TRIQS_PATH}/share/triqs/cmake) +# Append Pomerol installed files to the cmake load path +list(APPEND CMAKE_MODULE_PATH ${POMEROL_PATH}/share/pomerol/cmake) + +# Load TRIQS, including all predefined variables from TRIQS installation +find_package(TRIQS REQUIRED) + +if(NOT TRIQS_VERSION EQUAL 1.5) + message(FATAL_ERROR "The application requires the TRIQS library version 1.5 (got ${TRIQS_VERSION})") +endif() + +# Load Pomerol +find_package(pomerol REQUIRED) + +# Find MPI +find_package(MPI) + +# Get hash +triqs_get_git_hash(${CMAKE_SOURCE_DIR} "POMEROL2TRIQS") +if(${GIT_RESULT} EQUAL 0) + message(STATUS "Hash: ${POMEROL2TRIQS_GIT_HASH}") +endif(${GIT_RESULT} EQUAL 0) + +# We want to be installed in the TRIQS tree +set(CMAKE_INSTALL_PREFIX ${TRIQS_PATH}) + +message(STATUS "TRIQS : Adding compilation flags detected by the library") +add_definitions(${TRIQS_CXX_DEFINITIONS}) + +option(Tests "Enable Tests" ON) + +include_directories(c++) + +# Compile C++ code +add_subdirectory(c++) + +# Python interface +if (${TRIQS_WITH_PYTHON_SUPPORT}) + add_subdirectory(python) + if (${Tests}) + enable_testing() + add_subdirectory(test) + endif() +endif() diff --git a/LICENSE b/LICENSE new file mode 100644 index 0000000..94a9ed0 --- /dev/null +++ b/LICENSE @@ -0,0 +1,674 @@ + GNU GENERAL PUBLIC LICENSE + Version 3, 29 June 2007 + + Copyright (C) 2007 Free Software Foundation, Inc. + Everyone is permitted to copy and distribute verbatim copies + of this license document, but changing it is not allowed. + + Preamble + + The GNU General Public License is a free, copyleft license for +software and other kinds of works. + + The licenses for most software and other practical works are designed +to take away your freedom to share and change the works. By contrast, +the GNU General Public License is intended to guarantee your freedom to +share and change all versions of a program--to make sure it remains free +software for all its users. We, the Free Software Foundation, use the +GNU General Public License for most of our software; it applies also to +any other work released this way by its authors. You can apply it to +your programs, too. + + When we speak of free software, we are referring to freedom, not +price. Our General Public Licenses are designed to make sure that you +have the freedom to distribute copies of free software (and charge for +them if you wish), that you receive source code or can get it if you +want it, that you can change the software or use pieces of it in new +free programs, and that you know you can do these things. + + To protect your rights, we need to prevent others from denying you +these rights or asking you to surrender the rights. Therefore, you have +certain responsibilities if you distribute copies of the software, or if +you modify it: responsibilities to respect the freedom of others. + + For example, if you distribute copies of such a program, whether +gratis or for a fee, you must pass on to the recipients the same +freedoms that you received. You must make sure that they, too, receive +or can get the source code. And you must show them these terms so they +know their rights. + + Developers that use the GNU GPL protect your rights with two steps: +(1) assert copyright on the software, and (2) offer you this License +giving you legal permission to copy, distribute and/or modify it. + + For the developers' and authors' protection, the GPL clearly explains +that there is no warranty for this free software. For both users' and +authors' sake, the GPL requires that modified versions be marked as +changed, so that their problems will not be attributed erroneously to +authors of previous versions. + + Some devices are designed to deny users access to install or run +modified versions of the software inside them, although the manufacturer +can do so. This is fundamentally incompatible with the aim of +protecting users' freedom to change the software. The systematic +pattern of such abuse occurs in the area of products for individuals to +use, which is precisely where it is most unacceptable. Therefore, we +have designed this version of the GPL to prohibit the practice for those +products. If such problems arise substantially in other domains, we +stand ready to extend this provision to those domains in future versions +of the GPL, as needed to protect the freedom of users. + + Finally, every program is threatened constantly by software patents. +States should not allow patents to restrict development and use of +software on general-purpose computers, but in those that do, we wish to +avoid the special danger that patents applied to a free program could +make it effectively proprietary. To prevent this, the GPL assures that +patents cannot be used to render the program non-free. + + The precise terms and conditions for copying, distribution and +modification follow. + + TERMS AND CONDITIONS + + 0. Definitions. + + "This License" refers to version 3 of the GNU General Public License. + + "Copyright" also means copyright-like laws that apply to other kinds of +works, such as semiconductor masks. + + "The Program" refers to any copyrightable work licensed under this +License. Each licensee is addressed as "you". "Licensees" and +"recipients" may be individuals or organizations. + + To "modify" a work means to copy from or adapt all or part of the work +in a fashion requiring copyright permission, other than the making of an +exact copy. The resulting work is called a "modified version" of the +earlier work or a work "based on" the earlier work. + + A "covered work" means either the unmodified Program or a work based +on the Program. + + To "propagate" a work means to do anything with it that, without +permission, would make you directly or secondarily liable for +infringement under applicable copyright law, except executing it on a +computer or modifying a private copy. Propagation includes copying, +distribution (with or without modification), making available to the +public, and in some countries other activities as well. + + To "convey" a work means any kind of propagation that enables other +parties to make or receive copies. Mere interaction with a user through +a computer network, with no transfer of a copy, is not conveying. + + An interactive user interface displays "Appropriate Legal Notices" +to the extent that it includes a convenient and prominently visible +feature that (1) displays an appropriate copyright notice, and (2) +tells the user that there is no warranty for the work (except to the +extent that warranties are provided), that licensees may convey the +work under this License, and how to view a copy of this License. If +the interface presents a list of user commands or options, such as a +menu, a prominent item in the list meets this criterion. + + 1. Source Code. + + The "source code" for a work means the preferred form of the work +for making modifications to it. "Object code" means any non-source +form of a work. + + A "Standard Interface" means an interface that either is an official +standard defined by a recognized standards body, or, in the case of +interfaces specified for a particular programming language, one that +is widely used among developers working in that language. + + The "System Libraries" of an executable work include anything, other +than the work as a whole, that (a) is included in the normal form of +packaging a Major Component, but which is not part of that Major +Component, and (b) serves only to enable use of the work with that +Major Component, or to implement a Standard Interface for which an +implementation is available to the public in source code form. A +"Major Component", in this context, means a major essential component +(kernel, window system, and so on) of the specific operating system +(if any) on which the executable work runs, or a compiler used to +produce the work, or an object code interpreter used to run it. + + The "Corresponding Source" for a work in object code form means all +the source code needed to generate, install, and (for an executable +work) run the object code and to modify the work, including scripts to +control those activities. However, it does not include the work's +System Libraries, or general-purpose tools or generally available free +programs which are used unmodified in performing those activities but +which are not part of the work. For example, Corresponding Source +includes interface definition files associated with source files for +the work, and the source code for shared libraries and dynamically +linked subprograms that the work is specifically designed to require, +such as by intimate data communication or control flow between those +subprograms and other parts of the work. + + The Corresponding Source need not include anything that users +can regenerate automatically from other parts of the Corresponding +Source. + + The Corresponding Source for a work in source code form is that +same work. + + 2. Basic Permissions. + + All rights granted under this License are granted for the term of +copyright on the Program, and are irrevocable provided the stated +conditions are met. This License explicitly affirms your unlimited +permission to run the unmodified Program. The output from running a +covered work is covered by this License only if the output, given its +content, constitutes a covered work. This License acknowledges your +rights of fair use or other equivalent, as provided by copyright law. + + You may make, run and propagate covered works that you do not +convey, without conditions so long as your license otherwise remains +in force. You may convey covered works to others for the sole purpose +of having them make modifications exclusively for you, or provide you +with facilities for running those works, provided that you comply with +the terms of this License in conveying all material for which you do +not control copyright. Those thus making or running the covered works +for you must do so exclusively on your behalf, under your direction +and control, on terms that prohibit them from making any copies of +your copyrighted material outside their relationship with you. + + Conveying under any other circumstances is permitted solely under +the conditions stated below. Sublicensing is not allowed; section 10 +makes it unnecessary. + + 3. Protecting Users' Legal Rights From Anti-Circumvention Law. + + No covered work shall be deemed part of an effective technological +measure under any applicable law fulfilling obligations under article +11 of the WIPO copyright treaty adopted on 20 December 1996, or +similar laws prohibiting or restricting circumvention of such +measures. + + When you convey a covered work, you waive any legal power to forbid +circumvention of technological measures to the extent such circumvention +is effected by exercising rights under this License with respect to +the covered work, and you disclaim any intention to limit operation or +modification of the work as a means of enforcing, against the work's +users, your or third parties' legal rights to forbid circumvention of +technological measures. + + 4. Conveying Verbatim Copies. + + You may convey verbatim copies of the Program's source code as you +receive it, in any medium, provided that you conspicuously and +appropriately publish on each copy an appropriate copyright notice; +keep intact all notices stating that this License and any +non-permissive terms added in accord with section 7 apply to the code; +keep intact all notices of the absence of any warranty; and give all +recipients a copy of this License along with the Program. + + You may charge any price or no price for each copy that you convey, +and you may offer support or warranty protection for a fee. + + 5. Conveying Modified Source Versions. + + You may convey a work based on the Program, or the modifications to +produce it from the Program, in the form of source code under the +terms of section 4, provided that you also meet all of these conditions: + + a) The work must carry prominent notices stating that you modified + it, and giving a relevant date. + + b) The work must carry prominent notices stating that it is + released under this License and any conditions added under section + 7. This requirement modifies the requirement in section 4 to + "keep intact all notices". + + c) You must license the entire work, as a whole, under this + License to anyone who comes into possession of a copy. This + License will therefore apply, along with any applicable section 7 + additional terms, to the whole of the work, and all its parts, + regardless of how they are packaged. This License gives no + permission to license the work in any other way, but it does not + invalidate such permission if you have separately received it. + + d) If the work has interactive user interfaces, each must display + Appropriate Legal Notices; however, if the Program has interactive + interfaces that do not display Appropriate Legal Notices, your + work need not make them do so. + + A compilation of a covered work with other separate and independent +works, which are not by their nature extensions of the covered work, +and which are not combined with it such as to form a larger program, +in or on a volume of a storage or distribution medium, is called an +"aggregate" if the compilation and its resulting copyright are not +used to limit the access or legal rights of the compilation's users +beyond what the individual works permit. Inclusion of a covered work +in an aggregate does not cause this License to apply to the other +parts of the aggregate. + + 6. Conveying Non-Source Forms. + + You may convey a covered work in object code form under the terms +of sections 4 and 5, provided that you also convey the +machine-readable Corresponding Source under the terms of this License, +in one of these ways: + + a) Convey the object code in, or embodied in, a physical product + (including a physical distribution medium), accompanied by the + Corresponding Source fixed on a durable physical medium + customarily used for software interchange. + + b) Convey the object code in, or embodied in, a physical product + (including a physical distribution medium), accompanied by a + written offer, valid for at least three years and valid for as + long as you offer spare parts or customer support for that product + model, to give anyone who possesses the object code either (1) a + copy of the Corresponding Source for all the software in the + product that is covered by this License, on a durable physical + medium customarily used for software interchange, for a price no + more than your reasonable cost of physically performing this + conveying of source, or (2) access to copy the + Corresponding Source from a network server at no charge. + + c) Convey individual copies of the object code with a copy of the + written offer to provide the Corresponding Source. This + alternative is allowed only occasionally and noncommercially, and + only if you received the object code with such an offer, in accord + with subsection 6b. + + d) Convey the object code by offering access from a designated + place (gratis or for a charge), and offer equivalent access to the + Corresponding Source in the same way through the same place at no + further charge. You need not require recipients to copy the + Corresponding Source along with the object code. If the place to + copy the object code is a network server, the Corresponding Source + may be on a different server (operated by you or a third party) + that supports equivalent copying facilities, provided you maintain + clear directions next to the object code saying where to find the + Corresponding Source. Regardless of what server hosts the + Corresponding Source, you remain obligated to ensure that it is + available for as long as needed to satisfy these requirements. + + e) Convey the object code using peer-to-peer transmission, provided + you inform other peers where the object code and Corresponding + Source of the work are being offered to the general public at no + charge under subsection 6d. + + A separable portion of the object code, whose source code is excluded +from the Corresponding Source as a System Library, need not be +included in conveying the object code work. + + A "User Product" is either (1) a "consumer product", which means any +tangible personal property which is normally used for personal, family, +or household purposes, or (2) anything designed or sold for incorporation +into a dwelling. In determining whether a product is a consumer product, +doubtful cases shall be resolved in favor of coverage. For a particular +product received by a particular user, "normally used" refers to a +typical or common use of that class of product, regardless of the status +of the particular user or of the way in which the particular user +actually uses, or expects or is expected to use, the product. A product +is a consumer product regardless of whether the product has substantial +commercial, industrial or non-consumer uses, unless such uses represent +the only significant mode of use of the product. + + "Installation Information" for a User Product means any methods, +procedures, authorization keys, or other information required to install +and execute modified versions of a covered work in that User Product from +a modified version of its Corresponding Source. The information must +suffice to ensure that the continued functioning of the modified object +code is in no case prevented or interfered with solely because +modification has been made. + + If you convey an object code work under this section in, or with, or +specifically for use in, a User Product, and the conveying occurs as +part of a transaction in which the right of possession and use of the +User Product is transferred to the recipient in perpetuity or for a +fixed term (regardless of how the transaction is characterized), the +Corresponding Source conveyed under this section must be accompanied +by the Installation Information. But this requirement does not apply +if neither you nor any third party retains the ability to install +modified object code on the User Product (for example, the work has +been installed in ROM). + + The requirement to provide Installation Information does not include a +requirement to continue to provide support service, warranty, or updates +for a work that has been modified or installed by the recipient, or for +the User Product in which it has been modified or installed. Access to a +network may be denied when the modification itself materially and +adversely affects the operation of the network or violates the rules and +protocols for communication across the network. + + Corresponding Source conveyed, and Installation Information provided, +in accord with this section must be in a format that is publicly +documented (and with an implementation available to the public in +source code form), and must require no special password or key for +unpacking, reading or copying. + + 7. Additional Terms. + + "Additional permissions" are terms that supplement the terms of this +License by making exceptions from one or more of its conditions. +Additional permissions that are applicable to the entire Program shall +be treated as though they were included in this License, to the extent +that they are valid under applicable law. If additional permissions +apply only to part of the Program, that part may be used separately +under those permissions, but the entire Program remains governed by +this License without regard to the additional permissions. + + When you convey a copy of a covered work, you may at your option +remove any additional permissions from that copy, or from any part of +it. (Additional permissions may be written to require their own +removal in certain cases when you modify the work.) You may place +additional permissions on material, added by you to a covered work, +for which you have or can give appropriate copyright permission. + + Notwithstanding any other provision of this License, for material you +add to a covered work, you may (if authorized by the copyright holders of +that material) supplement the terms of this License with terms: + + a) Disclaiming warranty or limiting liability differently from the + terms of sections 15 and 16 of this License; or + + b) Requiring preservation of specified reasonable legal notices or + author attributions in that material or in the Appropriate Legal + Notices displayed by works containing it; or + + c) Prohibiting misrepresentation of the origin of that material, or + requiring that modified versions of such material be marked in + reasonable ways as different from the original version; or + + d) Limiting the use for publicity purposes of names of licensors or + authors of the material; or + + e) Declining to grant rights under trademark law for use of some + trade names, trademarks, or service marks; or + + f) Requiring indemnification of licensors and authors of that + material by anyone who conveys the material (or modified versions of + it) with contractual assumptions of liability to the recipient, for + any liability that these contractual assumptions directly impose on + those licensors and authors. + + All other non-permissive additional terms are considered "further +restrictions" within the meaning of section 10. If the Program as you +received it, or any part of it, contains a notice stating that it is +governed by this License along with a term that is a further +restriction, you may remove that term. If a license document contains +a further restriction but permits relicensing or conveying under this +License, you may add to a covered work material governed by the terms +of that license document, provided that the further restriction does +not survive such relicensing or conveying. + + If you add terms to a covered work in accord with this section, you +must place, in the relevant source files, a statement of the +additional terms that apply to those files, or a notice indicating +where to find the applicable terms. + + Additional terms, permissive or non-permissive, may be stated in the +form of a separately written license, or stated as exceptions; +the above requirements apply either way. + + 8. Termination. + + You may not propagate or modify a covered work except as expressly +provided under this License. Any attempt otherwise to propagate or +modify it is void, and will automatically terminate your rights under +this License (including any patent licenses granted under the third +paragraph of section 11). + + However, if you cease all violation of this License, then your +license from a particular copyright holder is reinstated (a) +provisionally, unless and until the copyright holder explicitly and +finally terminates your license, and (b) permanently, if the copyright +holder fails to notify you of the violation by some reasonable means +prior to 60 days after the cessation. + + Moreover, your license from a particular copyright holder is +reinstated permanently if the copyright holder notifies you of the +violation by some reasonable means, this is the first time you have +received notice of violation of this License (for any work) from that +copyright holder, and you cure the violation prior to 30 days after +your receipt of the notice. + + Termination of your rights under this section does not terminate the +licenses of parties who have received copies or rights from you under +this License. If your rights have been terminated and not permanently +reinstated, you do not qualify to receive new licenses for the same +material under section 10. + + 9. Acceptance Not Required for Having Copies. + + You are not required to accept this License in order to receive or +run a copy of the Program. Ancillary propagation of a covered work +occurring solely as a consequence of using peer-to-peer transmission +to receive a copy likewise does not require acceptance. However, +nothing other than this License grants you permission to propagate or +modify any covered work. These actions infringe copyright if you do +not accept this License. Therefore, by modifying or propagating a +covered work, you indicate your acceptance of this License to do so. + + 10. Automatic Licensing of Downstream Recipients. + + Each time you convey a covered work, the recipient automatically +receives a license from the original licensors, to run, modify and +propagate that work, subject to this License. You are not responsible +for enforcing compliance by third parties with this License. + + An "entity transaction" is a transaction transferring control of an +organization, or substantially all assets of one, or subdividing an +organization, or merging organizations. If propagation of a covered +work results from an entity transaction, each party to that +transaction who receives a copy of the work also receives whatever +licenses to the work the party's predecessor in interest had or could +give under the previous paragraph, plus a right to possession of the +Corresponding Source of the work from the predecessor in interest, if +the predecessor has it or can get it with reasonable efforts. + + You may not impose any further restrictions on the exercise of the +rights granted or affirmed under this License. For example, you may +not impose a license fee, royalty, or other charge for exercise of +rights granted under this License, and you may not initiate litigation +(including a cross-claim or counterclaim in a lawsuit) alleging that +any patent claim is infringed by making, using, selling, offering for +sale, or importing the Program or any portion of it. + + 11. Patents. + + A "contributor" is a copyright holder who authorizes use under this +License of the Program or a work on which the Program is based. The +work thus licensed is called the contributor's "contributor version". + + A contributor's "essential patent claims" are all patent claims +owned or controlled by the contributor, whether already acquired or +hereafter acquired, that would be infringed by some manner, permitted +by this License, of making, using, or selling its contributor version, +but do not include claims that would be infringed only as a +consequence of further modification of the contributor version. For +purposes of this definition, "control" includes the right to grant +patent sublicenses in a manner consistent with the requirements of +this License. + + Each contributor grants you a non-exclusive, worldwide, royalty-free +patent license under the contributor's essential patent claims, to +make, use, sell, offer for sale, import and otherwise run, modify and +propagate the contents of its contributor version. + + In the following three paragraphs, a "patent license" is any express +agreement or commitment, however denominated, not to enforce a patent +(such as an express permission to practice a patent or covenant not to +sue for patent infringement). To "grant" such a patent license to a +party means to make such an agreement or commitment not to enforce a +patent against the party. + + If you convey a covered work, knowingly relying on a patent license, +and the Corresponding Source of the work is not available for anyone +to copy, free of charge and under the terms of this License, through a +publicly available network server or other readily accessible means, +then you must either (1) cause the Corresponding Source to be so +available, or (2) arrange to deprive yourself of the benefit of the +patent license for this particular work, or (3) arrange, in a manner +consistent with the requirements of this License, to extend the patent +license to downstream recipients. "Knowingly relying" means you have +actual knowledge that, but for the patent license, your conveying the +covered work in a country, or your recipient's use of the covered work +in a country, would infringe one or more identifiable patents in that +country that you have reason to believe are valid. + + If, pursuant to or in connection with a single transaction or +arrangement, you convey, or propagate by procuring conveyance of, a +covered work, and grant a patent license to some of the parties +receiving the covered work authorizing them to use, propagate, modify +or convey a specific copy of the covered work, then the patent license +you grant is automatically extended to all recipients of the covered +work and works based on it. + + A patent license is "discriminatory" if it does not include within +the scope of its coverage, prohibits the exercise of, or is +conditioned on the non-exercise of one or more of the rights that are +specifically granted under this License. You may not convey a covered +work if you are a party to an arrangement with a third party that is +in the business of distributing software, under which you make payment +to the third party based on the extent of your activity of conveying +the work, and under which the third party grants, to any of the +parties who would receive the covered work from you, a discriminatory +patent license (a) in connection with copies of the covered work +conveyed by you (or copies made from those copies), or (b) primarily +for and in connection with specific products or compilations that +contain the covered work, unless you entered into that arrangement, +or that patent license was granted, prior to 28 March 2007. + + Nothing in this License shall be construed as excluding or limiting +any implied license or other defenses to infringement that may +otherwise be available to you under applicable patent law. + + 12. No Surrender of Others' Freedom. + + If conditions are imposed on you (whether by court order, agreement or +otherwise) that contradict the conditions of this License, they do not +excuse you from the conditions of this License. If you cannot convey a +covered work so as to satisfy simultaneously your obligations under this +License and any other pertinent obligations, then as a consequence you may +not convey it at all. For example, if you agree to terms that obligate you +to collect a royalty for further conveying from those to whom you convey +the Program, the only way you could satisfy both those terms and this +License would be to refrain entirely from conveying the Program. + + 13. Use with the GNU Affero General Public License. + + Notwithstanding any other provision of this License, you have +permission to link or combine any covered work with a work licensed +under version 3 of the GNU Affero General Public License into a single +combined work, and to convey the resulting work. The terms of this +License will continue to apply to the part which is the covered work, +but the special requirements of the GNU Affero General Public License, +section 13, concerning interaction through a network will apply to the +combination as such. + + 14. Revised Versions of this License. + + The Free Software Foundation may publish revised and/or new versions of +the GNU General Public License from time to time. Such new versions will +be similar in spirit to the present version, but may differ in detail to +address new problems or concerns. + + Each version is given a distinguishing version number. If the +Program specifies that a certain numbered version of the GNU General +Public License "or any later version" applies to it, you have the +option of following the terms and conditions either of that numbered +version or of any later version published by the Free Software +Foundation. If the Program does not specify a version number of the +GNU General Public License, you may choose any version ever published +by the Free Software Foundation. + + If the Program specifies that a proxy can decide which future +versions of the GNU General Public License can be used, that proxy's +public statement of acceptance of a version permanently authorizes you +to choose that version for the Program. + + Later license versions may give you additional or different +permissions. However, no additional obligations are imposed on any +author or copyright holder as a result of your choosing to follow a +later version. + + 15. Disclaimer of Warranty. + + THERE IS NO WARRANTY FOR THE PROGRAM, TO THE EXTENT PERMITTED BY +APPLICABLE LAW. EXCEPT WHEN OTHERWISE STATED IN WRITING THE COPYRIGHT +HOLDERS AND/OR OTHER PARTIES PROVIDE THE PROGRAM "AS IS" WITHOUT WARRANTY +OF ANY KIND, EITHER EXPRESSED OR IMPLIED, INCLUDING, BUT NOT LIMITED TO, +THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR +PURPOSE. THE ENTIRE RISK AS TO THE QUALITY AND PERFORMANCE OF THE PROGRAM +IS WITH YOU. SHOULD THE PROGRAM PROVE DEFECTIVE, YOU ASSUME THE COST OF +ALL NECESSARY SERVICING, REPAIR OR CORRECTION. + + 16. Limitation of Liability. + + IN NO EVENT UNLESS REQUIRED BY APPLICABLE LAW OR AGREED TO IN WRITING +WILL ANY COPYRIGHT HOLDER, OR ANY OTHER PARTY WHO MODIFIES AND/OR CONVEYS +THE PROGRAM AS PERMITTED ABOVE, BE LIABLE TO YOU FOR DAMAGES, INCLUDING ANY +GENERAL, SPECIAL, INCIDENTAL OR CONSEQUENTIAL DAMAGES ARISING OUT OF THE +USE OR INABILITY TO USE THE PROGRAM (INCLUDING BUT NOT LIMITED TO LOSS OF +DATA OR DATA BEING RENDERED INACCURATE OR LOSSES SUSTAINED BY YOU OR THIRD +PARTIES OR A FAILURE OF THE PROGRAM TO OPERATE WITH ANY OTHER PROGRAMS), +EVEN IF SUCH HOLDER OR OTHER PARTY HAS BEEN ADVISED OF THE POSSIBILITY OF +SUCH DAMAGES. + + 17. Interpretation of Sections 15 and 16. + + If the disclaimer of warranty and limitation of liability provided +above cannot be given local legal effect according to their terms, +reviewing courts shall apply local law that most closely approximates +an absolute waiver of all civil liability in connection with the +Program, unless a warranty or assumption of liability accompanies a +copy of the Program in return for a fee. + + END OF TERMS AND CONDITIONS + + How to Apply These Terms to Your New Programs + + If you develop a new program, and you want it to be of the greatest +possible use to the public, the best way to achieve this is to make it +free software which everyone can redistribute and change under these terms. + + To do so, attach the following notices to the program. It is safest +to attach them to the start of each source file to most effectively +state the exclusion of warranty; and each file should have at least +the "copyright" line and a pointer to where the full notice is found. + + + Copyright (C) + + 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 . + +Also add information on how to contact you by electronic and paper mail. + + If the program does terminal interaction, make it output a short +notice like this when it starts in an interactive mode: + + Copyright (C) + This program comes with ABSOLUTELY NO WARRANTY; for details type `show w'. + This is free software, and you are welcome to redistribute it + under certain conditions; type `show c' for details. + +The hypothetical commands `show w' and `show c' should show the appropriate +parts of the General Public License. Of course, your program's commands +might be different; for a GUI interface, you would use an "about box". + + You should also get your employer (if you work as a programmer) or school, +if any, to sign a "copyright disclaimer" for the program, if necessary. +For more information on this, and how to apply and follow the GNU GPL, see +. + + The GNU General Public License does not permit incorporating your program +into proprietary programs. If your program is a subroutine library, you +may consider it more useful to permit linking proprietary applications with +the library. If this is what you want to do, use the GNU Lesser General +Public License instead of this License. But first, please read +. diff --git a/README.md b/README.md new file mode 100644 index 0000000..66285f0 --- /dev/null +++ b/README.md @@ -0,0 +1,45 @@ +pomerol2triqs +============= + +Quick and dirty TRIQS wrapper around the Pomerol exact diagonalization library + +To learn how to use this wrapper, see `example` subdir in the source directory. + +Features +-------- + +* Diagonalization of finite fermionic models with Hamiltonians written in terms of second quantization operators. +* Calculation of single-particle Green's functions: `G(\tau)`, `G(i\omega_n)`, `G(\omega)`. +* Calculation of two-particle Green's functions: `G(\omega;\nu,\nu')` and `G(\omega;\ell,\ell')`. + +Notation for the two-particle Green's functions is adopted from the +[PhD thesis of Lewin Boehnke](http://ediss.sub.uni-hamburg.de/volltexte/2015/7325/pdf/Dissertation.pdf). + +Installation +------------ + +- Install the latest version of [Pomerol](http://aeantipov.github.io/pomerol/) exact diagonalization library (`master` branch). +- Install the latest version of [TRIQS](https://triqs.ipht.cnrs.fr/1.x/install.html) library (**`unstable`** branch). +- `git clone https://github.com/krivenko/pomerol2triqs.git pomerol2triqs.git` +- `mkdir pomerol2triqs.build && cd pomerol2triqs.build` +- `cmake ../pomerol2triqs.git -DCMAKE_BUILD_TYPE=Release -DTRIQS_PATH= -DPOMEROL_PATH=` +- `make` +- `make install` + +License +------- + +Copyright (C) 2017 Igor Krivenko + +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 . diff --git a/c++/CMakeLists.txt b/c++/CMakeLists.txt new file mode 100644 index 0000000..9c873d7 --- /dev/null +++ b/c++/CMakeLists.txt @@ -0,0 +1,8 @@ +add_library(pomerol2triqs_c pomerol_ed.cpp) + +target_link_libraries(pomerol2triqs_c ${TRIQS_LIBRARY_ALL} ${pomerol_LIBRARIES}) +include_directories(${TRIQS_INCLUDE_ALL} ${pomerol_INCLUDE_DIRS}) +triqs_set_rpath_for_target(pomerol2triqs_c) + +install(TARGETS pomerol2triqs_c DESTINATION lib) +install(FILES pomerol_ed.hpp DESTINATION include/pomerol2triqs) diff --git a/c++/g2_parameters.hpp b/c++/g2_parameters.hpp new file mode 100644 index 0000000..c828371 --- /dev/null +++ b/c++/g2_parameters.hpp @@ -0,0 +1,74 @@ +#pragma once + +#include +using triqs::hilbert_space::gf_struct_t; + +namespace pomerol2triqs { + + enum block_order_t { AABB, ABBA }; + enum channel_t { PP, PH, AllFermionic }; + + using g2_blocks_t = std::set>; + + struct g2_iw_inu_inup_params_t { + + /// Structure of G^2 blocks. + gf_struct_t gf_struct; + + /// Inverse temperature + double beta; + + /// Channel in which Matsubara frequency representation is defined. + channel_t channel = PH; + + /// Order of block indices in the definition of G^2. + block_order_t block_order = AABB; + + /// List of block index pairs of G^2 to measure. + /// default: measure all blocks + g2_blocks_t blocks = g2_blocks_t{}; + + /// Number of bosonic Matsubara frequencies. + int n_iw = 30; + + /// Number of fermionic Matsubara frequencies. + int n_inu = 30; + + g2_iw_inu_inup_params_t() {} + g2_iw_inu_inup_params_t(gf_struct_t const &gf_struct, double beta) : gf_struct(gf_struct), beta(beta) {} + }; + + struct g2_iw_l_lp_params_t { + + /// Structure of G^2 blocks. + gf_struct_t gf_struct; + + /// Inverse temperature + double beta; + + /// Channel in which Matsubara frequency representation is defined. + channel_t channel = PH; + + /// Order of block indices in the definition of G^2. + block_order_t block_order = AABB; + + /// List of block index pairs of G^2 to measure. + /// default: measure all blocks + g2_blocks_t blocks = g2_blocks_t{}; + + /// Number of bosonic Matsubara frequencies. + int n_iw = 30; + + /// Number of Legendre coefficients. + int n_l = 20; + + /// Maximum number of positive Matsubara frequencies in summation. + int n_inu_sum = 500; + + /// Tolerance for Matsubara frequency summation. + double inu_sum_tol = 1e-6; + + g2_iw_l_lp_params_t() {} + g2_iw_l_lp_params_t(gf_struct_t const &gf_struct, double beta) : gf_struct(gf_struct), beta(beta) {} + }; +} diff --git a/c++/pomerol_ed.cpp b/c++/pomerol_ed.cpp new file mode 100644 index 0000000..d6be072 --- /dev/null +++ b/c++/pomerol_ed.cpp @@ -0,0 +1,538 @@ +#include "pomerol_ed.hpp" + +#include +#include +#include + +namespace pomerol2triqs { + + // Generalization of triqs::utility::legendre_T + // See eq. (4.5) in Lewin Boehnke's thesis + inline std::complex t_bar(int o, int l) { + if (o == 0) return l == 0 ? 1 : 0; + + const double pi = boost::math::constants::pi(); + bool neg_o = false; + if (o < 0) { + neg_o = true; + o = -o; + } + + std::complex res = (sqrt(2 * l + 1) / sqrt(o)) * std::pow(1_j, o + l) * boost::math::cyl_bessel_j(l + 0.5, o * pi / 2); + // \bar T_{-ol} = \bar T_{ol}^* + return neg_o ? std::conj(res) : res; + } + + Pomerol::Lattice pomerol_ed::init() { + + Pomerol::Lattice l; + + std::map site_max_orb; + for (auto const &ind : index_converter) { + std::string pomerol_site; + int pomerol_orb; + std::tie(pomerol_site, pomerol_orb, std::ignore) = ind.second; + + auto it = site_max_orb.find(pomerol_site); + if (it == site_max_orb.end()) + site_max_orb[pomerol_site] = pomerol_orb; + else + it->second = std::max(it->second, pomerol_orb); + } + + for (auto const &site_orb : site_max_orb) l.addSite(new Pomerol::Lattice::Site(site_orb.first, site_orb.second + 1, 2)); + + return l; + } + + pomerol_ed::pomerol_ed(index_converter_t const &index_converter, bool verbose) + : verbose(verbose), index_converter(index_converter), bare_lattice(init()), index_info(bare_lattice.getSiteMap()) { + index_info.prepare(); + if (verbose && !comm.rank()) { + std::cout << "Pomerol: lattice sites" << std::endl; + bare_lattice.printSites(); + std::cout << "Pomerol: operator indices" << std::endl; + index_info.printIndices(); + } + } + + double pomerol_ed::diagonalize_prepare(many_body_op_t const &hamiltonian) { + // Workaround for the broken std::vector + struct bool_ { + bool b; + }; + + std::vector OperatorSequence; + std::vector SiteLabels; + std::vector Orbitals; + std::vector Spins; + + lattice.reset(new Pomerol::Lattice(bare_lattice)); + + double gs_shift = 0; + + for (auto const &term : hamiltonian) { + if (term.monomial.empty()) { + gs_shift = std::real(term.coef); + continue; // Constant term is unphysical anyway ... + } + OperatorSequence.clear(); + SiteLabels.clear(); + Orbitals.clear(); + Spins.clear(); + + for (auto o : term.monomial) { + auto it = index_converter.find(o.indices); + if (it == index_converter.end()) TRIQS_RUNTIME_ERROR << "diagonalize: invalid Hamiltonian, unexpected operator indices " << o.indices; + + OperatorSequence.push_back({o.dagger}); + + std::string site; + unsigned short orb; + Pomerol::spin s; + std::tie(site, orb, s) = it->second; + SiteLabels.push_back(site); + Orbitals.push_back(orb); + Spins.push_back(s); + } + + lattice->addTerm(new Pomerol::Lattice::Term(term.monomial.size(), reinterpret_cast(OperatorSequence.data()), term.coef, + SiteLabels.data(), Orbitals.data(), Spins.data())); + } + + storage.reset(new Pomerol::IndexHamiltonian(lattice.get(), index_info)); + storage->prepare(); + + if (verbose && !comm.rank()) { + std::cout << "Pomerol: terms of Hamiltonian" << std::endl; + std::cout << *storage << std::endl; + } + + return gs_shift; + } + + void pomerol_ed::diagonalize_main(double gs_shift) { + + // Classify many-body states + states_class.reset(new Pomerol::StatesClassification(index_info, *symm)); + states_class->compute(); + + // Matrix representation of the Hamiltonian + matrix_h.reset(new Pomerol::Hamiltonian(index_info, *storage, *states_class)); + matrix_h->prepare(comm); + matrix_h->compute(comm); + + // Get ground state energy + if (verbose && !comm.rank()) { std::cout << "Pomerol: ground state energy is " << matrix_h->getGroundEnergy() + gs_shift << std::endl; } + + // Reset containers, we will compute them later if needed + rho.release(); + ops_container.release(); + } + + void pomerol_ed::diagonalize(many_body_op_t const &hamiltonian, bool ignore_symmetries) { + + double gs_shift = diagonalize_prepare(hamiltonian); + + // Check the Hamiltonian commutes with the total number of particles + Pomerol::OperatorPresets::N N(index_info.getIndexSize()); + if (!storage->commutes(N)) TRIQS_RUNTIME_ERROR << "diagonalize: Hamiltonian does not conserve the total number of particles"; + + // Construct Symmetrizer + symm.reset(new Pomerol::Symmetrizer(index_info, *storage)); + symm->compute(ignore_symmetries); + + diagonalize_main(gs_shift); + } + + void pomerol_ed::diagonalize(many_body_op_t const &hamiltonian, std::vector const& integrals_of_motion) { + + double gs_shift = diagonalize_prepare(hamiltonian); + + std::vector iom; + for(auto const& op : integrals_of_motion) { + Pomerol::Operator pom_op; + for (auto const &term : op) { + Pomerol::Operator pom_term; + pom_term += term.coef; + for (auto o : term.monomial) { + Pomerol::ParticleIndex pom_ind = lookup_pomerol_index(o.indices); + + if (pom_ind == -1) + TRIQS_RUNTIME_ERROR << "diagonalize: invalid integral of motion, unexpected operator indices " << o.indices; + + if(o.dagger) + pom_term *= Pomerol::OperatorPresets::c_dag(pom_ind); + else + pom_term *= Pomerol::OperatorPresets::c(pom_ind); + } + + pom_op += pom_term; + } + iom.push_back(pom_op); + } + + // Construct Symmetrizer + symm.reset(new Pomerol::Symmetrizer(index_info, *storage)); + symm->compute(iom); + + diagonalize_main(gs_shift); + } + + Pomerol::ParticleIndex pomerol_ed::lookup_pomerol_index(indices_t const &i) const { + auto it = index_converter.find(i); + if (it == index_converter.end()) return -1; + std::string site; + unsigned short orb; + Pomerol::spin s; + std::tie(site, orb, s) = it->second; + return index_info.getIndex(site, orb, s); + } + + // Translate gf_struct into a set of ParticleIndex + std::set pomerol_ed::gf_struct_to_pomerol_indices(gf_struct_t const &gf_struct) const { + std::set indices; + for (auto const &b : gf_struct) { + for (auto const &i : b.second) { + auto pom_ind = lookup_pomerol_index({b.first, i}); + if (pom_ind == -1) TRIQS_RUNTIME_ERROR << "gf_struct_to_pomerol_indices: unexpected GF index " << b.first << "," << i; + indices.insert(pom_ind); + } + } + return indices; + } + + // Create the Density Matrix. + void pomerol_ed::compute_rho(double beta) { + if (!states_class || !matrix_h) TRIQS_RUNTIME_ERROR << "compute_rho: internal error!"; + + if (!rho || rho->beta != beta) { + if (verbose && !comm.rank()) std::cout << "Pomerol: computing density matrix for \\beta = " << beta << std::endl; + rho.reset(new Pomerol::DensityMatrix(*states_class, *matrix_h, beta)); + rho->prepare(); + rho->compute(); + } + } + + void pomerol_ed::compute_field_operators(gf_struct_t const &gf_struct) { + if (!states_class || !matrix_h) TRIQS_RUNTIME_ERROR << "compute_field_operators: internal error!"; + + auto new_ops = gf_struct_to_pomerol_indices(gf_struct); + if (!ops_container || computed_ops != new_ops) { + if (verbose && !comm.rank()) { + std::cout << "Pomerol: computing field operators with indices "; + bool comma = false; + for (auto i : new_ops) { + std::cout << (comma ? ", " : "") << i; + comma = true; + } + std::cout << std::endl; + } + + ops_container.reset(new Pomerol::FieldOperatorContainer(index_info, *states_class, *matrix_h)); + ops_container->prepareAll(new_ops); + ops_container->computeAll(); + + computed_ops = new_ops; + } + } + + ////////////////////// + // Green's function // + ////////////////////// + + template + block_gf pomerol_ed::compute_gf(gf_struct_t const &gf_struct, gf_mesh const &mesh, Filler filler) const { + + if (!states_class || !matrix_h || !rho || !ops_container) TRIQS_RUNTIME_ERROR << "compute_gf: internal error!"; + + struct index_visitor { + std::vector indices; + void operator()(int i) { indices.push_back(std::to_string(i)); } + void operator()(std::string s) { indices.push_back(s); } + }; + + std::vector block_names; + std::vector> g_blocks; + + for (auto const &bl : gf_struct) { + block_names.push_back(bl.first); + int n = bl.second.size(); + + index_visitor iv; + for (auto &ind : bl.second) { apply_visitor(iv, ind); } + std::vector> indices{{iv.indices, iv.indices}}; + + g_blocks.push_back(gf{mesh, {n, n}, indices}); + auto &g = g_blocks.back(); + + for (int i1 : range(n)) { + Pomerol::ParticleIndex pom_i1 = lookup_pomerol_index({bl.first, bl.second[i1]}); + for (int i2 : range(n)) { + Pomerol::ParticleIndex pom_i2 = lookup_pomerol_index({bl.first, bl.second[i2]}); + + if (verbose && !comm.rank()) + std::cout << "fill_gf: Filling GF component (" << bl.first << "," << bl.second[i1] << ")(" << bl.first << "," << bl.second[i2] << ")" + << std::endl; + auto g_el = slice_target_to_scalar(g, i1, i2); + + Pomerol::GreensFunction pom_g(*states_class, *matrix_h, ops_container->getAnnihilationOperator(pom_i1), + ops_container->getCreationOperator(pom_i2), *rho); + pom_g.prepare(); + pom_g.compute(); + + filler(g_el, pom_g); + g_el.singularity()(1) = (i1 == i2) ? 1 : 0; + } + } + } + return make_block_gf(block_names, std::move(g_blocks)); + } + + block_gf pomerol_ed::G_iw(gf_struct_t const &gf_struct, double beta, int n_iw) { + if (!matrix_h) TRIQS_RUNTIME_ERROR << "G_iw: no Hamiltonian has been diagonalized"; + compute_rho(beta); + compute_field_operators(gf_struct); + + auto filler = [](gf_view g_el, Pomerol::GreensFunction const &pom_g) { + for (auto iw : g_el.mesh()) g_el[iw] = pom_g(std::complex(iw)); + }; + return compute_gf(gf_struct, {beta, Fermion, n_iw}, filler); + } + + block_gf pomerol_ed::G_tau(gf_struct_t const &gf_struct, double beta, int n_tau) { + if (!matrix_h) TRIQS_RUNTIME_ERROR << "G_tau: no Hamiltonian has been diagonalized"; + compute_rho(beta); + compute_field_operators(gf_struct); + + auto filler = [](gf_view g_el, Pomerol::GreensFunction const &pom_g) { + for (auto tau : g_el.mesh()) g_el[tau] = pom_g.of_tau(tau); + }; + return compute_gf(gf_struct, {beta, Fermion, n_tau}, filler); + } + + block_gf pomerol_ed::G_w(gf_struct_t const &gf_struct, double beta, std::pair const &energy_window, int n_w, + double im_shift) { + if (!matrix_h) TRIQS_RUNTIME_ERROR << "G_w: no Hamiltonian has been diagonalized"; + compute_rho(beta); + compute_field_operators(gf_struct); + + auto filler = [im_shift](gf_view g_el, Pomerol::GreensFunction const &pom_g) { + for (auto w : g_el.mesh()) g_el[w] = pom_g(double(w) + 1_j * im_shift); + }; + return compute_gf(gf_struct, {energy_window.first, energy_window.second, n_w}, filler); + } + + /////////////////////////////////// + // Two-particle Green's function // + /////////////////////////////////// + + template + auto pomerol_ed::compute_g2(gf_struct_t const &gf_struct, gf_mesh const &mesh, block_order_t block_order, g2_blocks_t const &g2_blocks, + Filler filler) const -> block2_gf> { + + if (!states_class || !matrix_h || !rho || !ops_container) TRIQS_RUNTIME_ERROR << "compute_g2: internal error!"; + + bool compute_all_blocks = g2_blocks.empty(); + + std::vector>>> gf_vecvec; + std::vector block_names; + + for (auto const &bl1 : gf_struct) { + auto &A = bl1.first; + int A_size = bl1.second.size(); + int s1 = A_size; + block_names.push_back(A); + + std::vector>> gf_vec; + for (auto const &bl2 : gf_struct) { + auto &B = bl2.first; + int B_size = bl2.second.size(); + int s3 = B_size; + + int s2 = block_order == AABB ? s1 : s3; + int s4 = block_order == AABB ? s3 : s1; + + gf_vec.emplace_back(mesh, make_shape(s1, s2, s3, s4)); + + if (compute_all_blocks || g2_blocks.count({A, B})) { + auto const &A_inner = bl1.second; + auto const &B_inner = bl2.second; + + auto &g2_block = gf_vec.back(); + + for (int a : range(A_size)) + for (int b : range(A_size)) + for (int c : range(B_size)) + for (int d : range(B_size)) { + + if (verbose && !comm.rank()) { + std::cout << "compute_g2: filling G^2 element "; + if (block_order == AABB) { + std::cout << "(" << A << "," << a << ")"; + std::cout << "(" << A << "," << b << ")"; + std::cout << "(" << B << "," << c << ")"; + std::cout << "(" << B << "," << d << ")"; + } else { + std::cout << "(" << A << "," << a << ")"; + std::cout << "(" << B << "," << d << ")"; + std::cout << "(" << B << "," << c << ")"; + std::cout << "(" << A << "," << b << ")"; + } + std::cout << std::endl; + } + + auto g2_el = block_order == AABB ? slice_target_to_scalar(g2_block, a, b, c, d) : slice_target_to_scalar(g2_block, a, d, c, b); + + Pomerol::ParticleIndex pom_i1 = lookup_pomerol_index({A, A_inner[b]}); + Pomerol::ParticleIndex pom_i2 = lookup_pomerol_index({B, B_inner[d]}); + Pomerol::ParticleIndex pom_i3 = lookup_pomerol_index({A, A_inner[a]}); + Pomerol::ParticleIndex pom_i4 = lookup_pomerol_index({B, B_inner[c]}); + + Pomerol::TwoParticleGF pom_g2(*states_class, *matrix_h, ops_container->getAnnihilationOperator(pom_i1), + ops_container->getAnnihilationOperator(pom_i2), ops_container->getCreationOperator(pom_i3), + ops_container->getCreationOperator(pom_i4), *rho); + pom_g2.prepare(); + pom_g2.compute(false, {}, comm); + + filler(g2_el, pom_g2); + } + } + } + gf_vecvec.emplace_back(std::move(gf_vec)); + } + + return make_block2_gf(block_names, block_names, std::move(gf_vecvec)); + } + + auto pomerol_ed::G2_iw_inu_inup(g2_iw_inu_inup_params_t const &p) -> block2_gf> { + if (!matrix_h) TRIQS_RUNTIME_ERROR << "G2_iw_inu_inup: no Hamiltonian has been diagonalized"; + compute_rho(p.beta); + compute_field_operators(p.gf_struct); + + if (verbose && !comm.rank()) std::cout << "G2_iw_inu_inup: filling output container" << std::endl; + + auto filler = [&p, this](gf_view g2_el, auto const &pom_g2) { + long mesh_index = 0; + for (auto w_nu_nup : g2_el.mesh()) { + if ((mesh_index++) % comm.size() != comm.rank()) continue; + + if (p.channel == AllFermionic) { + + int n1 = std::get<0>(w_nu_nup).index(); + int n2 = std::get<1>(w_nu_nup).index(); + int n3 = std::get<2>(w_nu_nup).index(); + + if (p.block_order == AABB) + g2_el[w_nu_nup] = -pom_g2(n2, n1 + n3 - n2, n1); + else + g2_el[w_nu_nup] = +pom_g2(n1 + n3 - n2, n2, n1); + + } else { // p.channel == PH or PP + + int w_n = std::get<0>(w_nu_nup).index(); + int nu_n = std::get<1>(w_nu_nup).index(); + int nup_n = std::get<2>(w_nu_nup).index(); + + int W_n = p.channel == PH ? w_n + nu_n : w_n - nup_n - 1; + + if (p.block_order == AABB) { + g2_el[w_nu_nup] = -pom_g2(W_n, nup_n, nu_n); + } else { + g2_el[w_nu_nup] = +pom_g2(nup_n, W_n, nu_n); + } + } + } + }; + + gf_mesh mesh_b{p.beta, Boson, p.n_iw}; + gf_mesh mesh_f{p.beta, Fermion, p.n_inu}; + + gf_mesh mesh_bff{mesh_b, mesh_f, mesh_f}; + gf_mesh mesh_fff{mesh_f, mesh_f, mesh_f}; + + block2_gf> g2; + + if (p.channel == AllFermionic) + g2 = compute_g2(p.gf_struct, mesh_fff, p.block_order, p.blocks, filler); + else + g2 = compute_g2(p.gf_struct, mesh_bff, p.block_order, p.blocks, filler); + + g2() = mpi_all_reduce(g2(), comm); + + return g2; + } + + auto pomerol_ed::G2_iw_l_lp(g2_iw_l_lp_params_t const &p) -> block2_gf> { + if (!matrix_h) TRIQS_RUNTIME_ERROR << "G2_iw_l_lp: no Hamiltonian has been diagonalized"; + compute_rho(p.beta); + compute_field_operators(p.gf_struct); + + gf_mesh mesh{{p.beta, Boson, p.n_iw}, {p.beta, Fermion, static_cast(p.n_l)}, {p.beta, Fermion, static_cast(p.n_l)}}; + + if (verbose && !comm.rank()) std::cout << "G2_iw_l_lp: filling output container" << std::endl; + + auto filler = [&p, this](gf_view g2_el, auto const &pom_g2) { + + auto get_g2_iw_inu_inup_val = [&p, &pom_g2](long w_m, long nu_n, long nup_n) { + int W_n = p.channel == PH ? w_m + nu_n : w_m - nup_n - 1; + if (p.block_order == AABB) + return -pom_g2(W_n, nup_n, nu_n); + else + return +pom_g2(nup_n, W_n, nu_n); + }; + + array, 2> border_contrib(p.n_l, p.n_l); + array llp_element_converged(p.n_l, p.n_l); + int n_llp_elements_converged; + + long mesh_index = 0; + for (auto iw : std::get<0>(g2_el.mesh())) { + if((mesh_index++) % comm.size() != comm.rank()) continue; + + int w_m = iw.index(); + + llp_element_converged() = false; + n_llp_elements_converged = 0; + + // Summation over n and n' is done by adding new border points to + // a square summation domain. + for (int r = 0; r < p.n_inu_sum; ++r) { // r is current size of the domain + if (n_llp_elements_converged == p.n_l * p.n_l) break; + + border_contrib() = 0; + for (int n = -r; n <= r; ++n) { + auto iw_val_1 = get_g2_iw_inu_inup_val(w_m, r, n); + auto iw_val_2 = get_g2_iw_inu_inup_val(w_m, -r - 1, n - 1); + auto iw_val_3 = get_g2_iw_inu_inup_val(w_m, n - 1, r); + auto iw_val_4 = get_g2_iw_inu_inup_val(w_m, n, -r - 1); + + for (auto l : std::get<1>(g2_el.mesh())) + for (auto lp : std::get<2>(g2_el.mesh())) { + if (llp_element_converged(l, lp)) continue; + + using std::conj; + std::complex val = 0; + val += t_bar(2 * r + w_m + 1, l) * iw_val_1 * conj(t_bar(2 * n + w_m + 1, lp)); + val += t_bar(2 * (-r - 1) + w_m + 1, l) * iw_val_2 * conj(t_bar(2 * (n - 1) + w_m + 1, lp)); + val += t_bar(2 * (n - 1) + w_m + 1, l) * iw_val_3 * conj(t_bar(2 * r + w_m + 1, lp)); + val += t_bar(2 * n + w_m + 1, l) * iw_val_4 * conj(t_bar(2 * (-r - 1) + w_m + 1, lp)); + + g2_el[{iw, l, lp}] += val; + border_contrib(l, lp) += val; + if (std::abs(border_contrib(l, lp)) < p.inu_sum_tol) { + llp_element_converged(l, lp) = true; + ++n_llp_elements_converged; + } + } + } + } + } + }; + + auto g2 = compute_g2(p.gf_struct, mesh, p.block_order, p.blocks, filler); + g2() = mpi_all_reduce(g2(), comm); + + return g2; + } +} diff --git a/c++/pomerol_ed.hpp b/c++/pomerol_ed.hpp new file mode 100644 index 0000000..1ca5633 --- /dev/null +++ b/c++/pomerol_ed.hpp @@ -0,0 +1,96 @@ +#pragma once + +#include +#include +#include +#include +#include + +#include +#include +#include +#include +#include +#include + +#include "g2_parameters.hpp" + +namespace pomerol2triqs { + +#ifdef POMEROL_COMPLEX_MATRIX_ELEMENTS + using h_scalar_t = std::complex; +#else + using h_scalar_t = double; +#endif + + // Operator with real or complex value + using many_body_op_t = triqs::operators::many_body_operator_generic; + + using namespace triqs::gfs; + using triqs::hilbert_space::gf_struct_t; + + using indices_t = triqs::hilbert_space::fundamental_operator_set::indices_t; + using pomerol_indices_t = std::tuple; + using index_converter_t = std::map; + + class pomerol_ed { + + triqs::mpi::communicator comm; + + const bool verbose; + index_converter_t index_converter; + Pomerol::Lattice bare_lattice; + Pomerol::IndexClassification index_info; + + std::unique_ptr lattice; + std::unique_ptr storage; + std::unique_ptr symm; + std::unique_ptr states_class; + std::unique_ptr matrix_h; + std::unique_ptr rho; + std::set computed_ops; + std::unique_ptr ops_container; + + Pomerol::Lattice init(); + Pomerol::ParticleIndex lookup_pomerol_index(indices_t const &i) const; + std::set gf_struct_to_pomerol_indices(gf_struct_t const &gf_struct) const; + double diagonalize_prepare(many_body_op_t const &hamiltonian); + void diagonalize_main(double gs_shift); + void compute_rho(double beta); + void compute_field_operators(gf_struct_t const &gf_struct); + template block_gf compute_gf(gf_struct_t const &gf_struct, gf_mesh const &mesh, Filler filler) const; + + using w_nu_nup_t = cartesian_product; + using w_l_lp_t = cartesian_product; + template + block2_gf> compute_g2(gf_struct_t const &gf_struct, gf_mesh const &mesh, block_order_t block_order, + g2_blocks_t const &g2_blocks, Filler filler) const; + + public: + /// Create a new solver object + pomerol_ed(index_converter_t const &index_converter, bool verbose = false); + + /// Diagonalize Hamiltonian optionally employing conservation of N and S_z + void diagonalize(many_body_op_t const &hamiltonian, bool ignore_symmetries = false); + + /// Diagonalize Hamiltonian using provided integrals of motion + void diagonalize(many_body_op_t const &hamiltonian, std::vector const& integrals_of_motion); + + /// Green's function in Matsubara frequencies + block_gf G_iw(gf_struct_t const &gf_struct, double beta, int n_iw); + + /// Green's function in imaginary time + block_gf G_tau(gf_struct_t const &gf_struct, double beta, int n_tau); + + /// Retarded Green's function on real energy axis + block_gf G_w(gf_struct_t const &gf_struct, double beta, std::pair const &energy_window, int n_w, double im_shift = 0); + + /// Two-particle Green's function, Matsubara frequencies + TRIQS_WRAP_ARG_AS_DICT + block2_gf> G2_iw_inu_inup(g2_iw_inu_inup_params_t const &p); + + /// Two-particle Green's function, bosonic Matsubara frequency + Legendre coefficients + TRIQS_WRAP_ARG_AS_DICT + block2_gf> G2_iw_l_lp(g2_iw_l_lp_params_t const &p); + }; +} diff --git a/cmake/sitecustomize.py b/cmake/sitecustomize.py new file mode 100644 index 0000000..0f31ba9 --- /dev/null +++ b/cmake/sitecustomize.py @@ -0,0 +1,8 @@ +def application_pytriqs_import(name,*args,**kwargs): + if name.startswith('@package_name@'): + name = name[len('@package_name@')+1:] + return builtin_import(name,*args,**kwargs) + +import __builtin__ +__builtin__.__import__, builtin_import = application_pytriqs_import, __builtin__.__import__ + diff --git a/doc/.gitignore b/doc/.gitignore new file mode 100644 index 0000000..6b7f2eb --- /dev/null +++ b/doc/.gitignore @@ -0,0 +1,6 @@ +*.aux +*.bbl +*.blg +*.pdf +*.out +*.synctex.gz diff --git a/doc/g2_conv.tex b/doc/g2_conv.tex new file mode 100644 index 0000000..de8637a --- /dev/null +++ b/doc/g2_conv.tex @@ -0,0 +1,134 @@ +\documentclass[a4paper,12pt]{article} +\usepackage[utf8]{inputenc} +\usepackage[margin=0.7in]{geometry} +\usepackage{amsmath,amssymb} +\usepackage{caption} +\usepackage{verbatim} +\usepackage{hyperref} +\DeclareMathAlphabet{\mathpzc}{OT1}{pzc}{m}{it} + +\newcommand{\aver}[1]{\ensuremath{\langle#1\rangle}} +\renewcommand{\t}{\ensuremath{\tau}} +\newcommand{\w}{\ensuremath{\omega}} +\newcommand{\W}{\ensuremath{\Omega}} +\newcommand{\n}{\ensuremath{\nu}} +\newcommand{\TT}{\ensuremath{\mathbb{T}_\t}} + +\newcommand{\pom}{\ensuremath{\mathtt{Pomerol}}} + +\begin{document} +\title{Relations between conventions for the two-particle Green's functions + used in \texttt{Pomerol} and \texttt{TRIQS/CTHYB}} +\author{Igor Krivenko} +\maketitle + +\texttt{Pomerol} definition of the two-particle Green's function: +\begin{equation}\label{g2_pomerol} + G^{(2),\pom}_{\alpha\beta\gamma\delta}(\w_1,\w_2,\w_3) = + \int_0^\beta d\t_1d\t_1d\t_3 + e^{i\w_1\t_1 + i\w_2\t_2 - i\w_3\t_3} + \aver{\TT + c_\alpha(\t_1)c_\beta(\t_2)c^\dag_\gamma(\t_3)c^\dag_\delta(0)}. +\end{equation} + +\texttt{TRIQS/CTHYB} definitions: +\begin{itemize} + \item \textit{All fermionic case}, + \begin{equation}\label{g2_allfermionic} + G^{(2),3\nu}_{\alpha\beta\gamma\delta}(\nu_1,\nu_2,\nu_3) = + \int_0^\beta d\t_1d\t_1d\t_3 + e^{-i\nu_1\t_1 + i\nu_2\t_2 - i\nu_3\t_3} + \aver{\TT + c^\dagger_\alpha(\t_1)c_\beta(\t_2)c^\dagger_\gamma(\t_3)c_\delta(0)}. + \end{equation} + \item \textit{Particle-hole channel}, + \begin{equation}\label{g2_ph} + G^{(2),ph}_{\alpha\beta\gamma\delta}(\w;\n,\n') = + \frac{1}{\beta}\int_0^\beta d\t_1d\t_2d\t_3d\t_4\ + e^{-i\n\t_1} e^{i(\n+\w)\tau_2} e^{-i(\n'+\w)\t_3} e^{i\n'\t_4} + \aver{\TT + c^\dag_\alpha(\t_1) c_\beta(\t_2) c^\dag_\gamma(\t_3) + c_\delta(\t_4)}. + \end{equation} + + \item \textit{Particle-particle channel}, + \begin{equation}\label{g2_pp} + G^{(2),pp}_{\alpha\beta\gamma\delta}(\w;\n,\n') = + \frac{1}{\beta}\int_0^\beta d\t_1d\t_2d\t_3d\t_4\ + e^{-i\n\t_1} e^{i(\w-\n')\t_2} e^{-i(\w-\n)\t_3} e^{i\n'\t_4} + \aver{\TT + c^\dag_\alpha(\t_1) c_\beta(\t_2) c^\dag_\gamma(\t_3) + c_\delta(\t_4)}. + \end{equation} +\end{itemize} + +Relation between ph- and pp-channels: +\begin{equation} + G^{(2),pp}_{\alpha\beta\gamma\delta}(\w;\n,\n') = + G^{(2),ph}_{\alpha\beta\gamma\delta}(\w-\n-\n';\n,\n'). +\end{equation} + +General \texttt{Pomerol}-\texttt{TRIQS/CTHYB} relations disregarding block +structure. +\begin{align} + G^{(2),3\nu}_{\alpha\beta\gamma\delta}(\nu_1,\nu_2,\nu_3) &= + -G^{(2),\pom}_{\beta\delta\alpha\gamma}(\nu_2,\nu_1+\nu_3-\nu_2,\nu_1),\\ + G^{(2),ph}_{\alpha\beta\gamma\delta}(\w;\n,\n') &= + -G^{(2),\pom}_{\beta\delta\alpha\gamma}(\w+\n,\n',\n),\\ + G^{(2),pp}_{\alpha\beta\gamma\delta}(\w;\n,\n') &= + -G^{(2),\pom}_{\beta\delta\alpha\gamma}(\w-\n',\n',\n). +\end{align} + +Relations between $AABB$ and $ABBA$ block structures of $G^{(2)}$ in +\texttt{TRIQS/CTHYB}. +\begin{align} +G^{(2),ph}_{(A,a)(B,d)(B,c)(A,b)}(\w;\n,\n') &= - +G^{(2),ph}_{(A,a)(A,b)(B,c)(B,d)}(\n'-\n;\n,\n+\w),\\ +G^{(2),pp}_{(A,a)(B,d)(B,c)(A,b)}(\w;\n,\n') &= - +G^{(2),pp}_{(A,a)(A,b)(B,c)(B,d)}(\w;\n,\w-\n'). +\end{align} + +Relations for $AABB$ block structure of $G^{(2)}$ in \texttt{TRIQS/CTHYB}. +\begin{align} + G^{(2),ph}_{(A,a)(A,b)(B,c)(B,d)}(\w;\n,\n') &= + -G^{(2),\pom}_{(A,b)(B,d)(A,a)(B,c)}(\w+\n,\n',\n),\\ + G^{(2),pp}_{(A,a)(A,b)(B,c)(B,d)}(\w;\n,\n') &= + -G^{(2),\pom}_{(A,b)(B,d)(A,a)(B,c)}(\w-\n',\n',\n). +\end{align} + +Relations for $ABBA$ block structure of $G^{(2)}$ in \texttt{TRIQS/CTHYB}. +\begin{align} + G^{(2),ph}_{(A,a)(B,d)(B,c)(A,b)}(\w;\n,\n') &= + G^{(2),\pom}_{(A,b)(B,d)(A,a)(B,c)}(\n',\w+\n,\n),\\ + G^{(2),pp}_{(A,a)(B,d)(B,c)(A,b)}(\w;\n,\n') &= + G^{(2),\pom}_{(A,b)(B,d)(A,a)(B,c)}(\n',\w-\n',\n). +\end{align} + +Mixed Matsubara/Legendre representation of $G^{(2)}$ as a transformation of +$G^{(2)}_{\alpha\beta\gamma\delta}(\w;\nu,\nu')$ (both channels) +\cite{LewinThesis}. +\begin{equation}\label{legendre_transform} + G^{(2)}_{\alpha\beta\gamma\delta}(\w;\ell,\ell') \equiv + \sum_{n,n'\in\mathbb{Z}} + \bar T_{2n+m+1,\ell} + G^{(2)}_{\alpha\beta\gamma\delta}(\w;\nu_n,\nu_{n'}) + \bar T^*_{2n'+m+1,\ell'}, +\end{equation} +\begin{equation} + \bar T_{o,\ell} \equiv \frac{\sqrt{2\ell+1}}{\beta} + \int_0^\beta d\t e^{io\pi\frac{\t}{\beta}} P_\ell[x(\t)] = + \sqrt{2\ell+1}i^o i^\ell j_\ell\left(\frac{o\pi}{2}\right). +\end{equation} + +\bibliographystyle{plain} + +\begin{thebibliography}{9} + \bibitem{LewinThesis} + Lewin Volker Boehnke, + \emph{Susceptibilities in materials with multiple strongly correlated + orbitals} + PhD thesis, Universit\"at Hamburg, 2015, + \url{http://ediss.sub.uni-hamburg.de/volltexte/2015/7325/pdf/Dissertation.pdf} +\end{thebibliography} + +\end{document} diff --git a/example/2band.atom.py b/example/2band.atom.py new file mode 100644 index 0000000..3840ddc --- /dev/null +++ b/example/2band.atom.py @@ -0,0 +1,159 @@ +from pytriqs.archive import HDFArchive +from pytriqs.gf import * +from pytriqs.operators import Operator, c, c_dag, n +from pytriqs.operators.util.hamiltonians import h_int_kanamori +from pytriqs.utility import mpi +from pytriqs.applications.impurity_solvers.pomerol2triqs import PomerolED +import numpy as np +from itertools import product + +# 2-orbital Hubbard-Kanamori atom + +#################### +# Input parameters # +#################### + +beta = 10.0 # Inverse temperature +num_orb = 2 # Number of orbitals +U = 2.0 # Coulomb repulsion +mu = 1.5 # Chemical potential +J = 0.2 # Hund coupling + +spin_names = ("up", "dn") +orb_names = range(num_orb) + +# Number of Matsubara frequencies for GF calculation +n_iw = 1024 +# Number of imaginary time slices for GF calculation +n_tau = 10001 + +# Energy window for real frequency GF calculation +energy_window = (-5, 5) +# Number of frequency points for real frequency GF calculation +n_w = 1000 + +# Number of bosonic Matsubara frequencies for G^2 calculations +g2_n_iw = 5 +# Number of fermionic Matsubara frequencies for G^2 calculations +g2_n_inu = 10 +# Number of Legendre coefficients for G^2 calculations +g2_n_l = 10 +# Block index combinations for G^2 calculations +g2_blocks = set([("up", "up"), ("up", "dn"), ("dn", "up")]) + +gf_struct = {"up" : orb_names, "dn" : orb_names} +print "Block structure of single-particle Green's functions:", gf_struct + +# Conversion from TRIQS to Pomerol notation for operator indices +# TRIQS: block_name, inner_index +# Pomerol: site_label, orbital_index, spin_name +index_converter = {(sn, o) : ("loc", o, "down" if sn == "dn" else "up") + for sn, o in product(spin_names, orb_names)} + +# Make PomerolED solver object +ed = PomerolED(index_converter, verbose = True) + +# Number of particles on the impurity +N = sum(n(sn, o) for sn, o in product(spin_names, orb_names)) + +# Hamiltonian +H = h_int_kanamori(spin_names, orb_names, + np.array([[0, U-3*J], [U-3*J, 0]]), + np.array([[U, U-2*J], [U-2*J, U]]), + J, True) +H -= mu*N + +# Diagonalize H +ed.diagonalize(H) + +# Compute G(i\omega) +G_iw = ed.G_iw(gf_struct, beta, n_iw) + +# Compute G(\tau) +G_tau = ed.G_tau(gf_struct, beta, n_tau) + +# Compute G(\omega) +G_w = ed.G_w(gf_struct, beta, energy_window, n_w, 0.01) + +########### +# G^{(2)} # +########### + +common_g2_params = {'gf_struct' : gf_struct, + 'beta' : beta, + 'blocks' : g2_blocks, + 'n_iw' : g2_n_iw} + +############################### +# G^{(2)}(i\omega;i\nu,i\nu') # +############################### + +# Compute G^{(2),ph}(i\omega;i\nu,i\nu'), AABB block order +G2_iw_inu_inup_ph_AABB = ed.G2_iw_inu_inup(channel = "PH", + block_order = "AABB", + n_inu = g2_n_inu, + **common_g2_params) + +# Compute G^{(2),ph}(i\omega;i\nu,i\nu'), ABBA block order +G2_iw_inu_inup_ph_ABBA = ed.G2_iw_inu_inup(channel = "PH", + block_order = "ABBA", + n_inu = g2_n_inu, + **common_g2_params) + +# Compute G^{(2),pp}(i\omega;i\nu,i\nu'), AABB block order +G2_iw_inu_inup_pp_AABB = ed.G2_iw_inu_inup(channel = "PP", + block_order = "AABB", + n_inu = g2_n_inu, + **common_g2_params) + +# Compute G^{(2),pp}(i\omega;i\nu,i\nu'), ABBA block order +G2_iw_inu_inup_pp_ABBA = ed.G2_iw_inu_inup(channel = "PP", + block_order = "ABBA", + n_inu = g2_n_inu, + **common_g2_params) + +######################### +# G^{(2)}(i\omega;l,l') # +######################### + +# Compute G^{(2),ph}(i\omega;l,l'), AABB block order +G2_iw_l_lp_ph_AABB = ed.G2_iw_l_lp(channel = "PH", + block_order = "AABB", + n_l = g2_n_l, + **common_g2_params) + +# Compute G^{(2),ph}(i\omega;l,l'), ABBA block order +G2_iw_l_lp_ph_ABBA = ed.G2_iw_l_lp(channel = "PH", + block_order = "ABBA", + n_l = g2_n_l, + **common_g2_params) + +# Compute G^{(2),pp}(i\omega;l,l'), AABB block order +G2_iw_l_lp_pp_AABB = ed.G2_iw_l_lp(channel = "PP", + block_order = "AABB", + n_l = g2_n_l, + **common_g2_params) + +# Compute G^{(2),pp}(i\omega;l,l'), ABBA block order +G2_iw_l_lp_pp_ABBA = ed.G2_iw_l_lp(channel = "PP", + block_order = "ABBA", + n_l = g2_n_l, + **common_g2_params) + +################ +# Save results # +################ + +if mpi.is_master_node(): + with HDFArchive('2band.atom.h5', 'w') as ar: + ar['G_iw'] = G_iw + ar['G_tau'] = G_tau + ar['G_w'] = G_w + ar['G2_iw_inu_inup_ph_AABB'] = G2_iw_inu_inup_ph_AABB + ar['G2_iw_inu_inup_ph_ABBA'] = G2_iw_inu_inup_ph_ABBA + ar['G2_iw_inu_inup_pp_AABB'] = G2_iw_inu_inup_pp_AABB + ar['G2_iw_inu_inup_pp_ABBA'] = G2_iw_inu_inup_pp_ABBA + ar['G2_iw_l_lp_ph_AABB'] = G2_iw_l_lp_ph_AABB + ar['G2_iw_l_lp_ph_ABBA'] = G2_iw_l_lp_ph_ABBA + ar['G2_iw_l_lp_pp_AABB'] = G2_iw_l_lp_pp_AABB + ar['G2_iw_l_lp_pp_ABBA'] = G2_iw_l_lp_pp_ABBA diff --git a/example/2band.py b/example/2band.py new file mode 100644 index 0000000..56d1f08 --- /dev/null +++ b/example/2band.py @@ -0,0 +1,182 @@ +from pytriqs.archive import HDFArchive +from pytriqs.gf import * +from pytriqs.operators import Operator, c, c_dag, n +from pytriqs.operators.util.hamiltonians import h_int_kanamori +from pytriqs.utility import mpi +from pytriqs.applications.impurity_solvers.pomerol2triqs import PomerolED +import numpy as np +from itertools import product + +# 2-orbital impurity Anderson model (bath: 1 site * 2 orbitals) + +#################### +# Input parameters # +#################### + +beta = 10.0 # Inverse temperature +num_orb = 2 # Number of orbitals +U = 2.0 # Coulomb repulsion +mu = 1.5 # Chemical potential +J = 0.2 # Hund coupling + +# Levels of the bath sites +epsilon = np.array([-0.2, 0.2]) +# Hopping matrix +V = 0.7*np.eye(num_orb) + 0.1*(np.ones((num_orb, num_orb)) - np.eye(num_orb)) + +spin_names = ("up", "dn") +orb_names = range(num_orb) + +# Number of Matsubara frequencies for GF calculation +n_iw = 1024 + +# Number of imaginary time slices for GF calculation +n_tau = 10001 + +# Energy window for real frequency GF calculation +energy_window = (-5, 5) +# Number of frequency points for real frequency GF calculation +n_w = 1000 + +# Number of bosonic Matsubara frequencies for G^2 calculations +g2_n_iw = 5 +# Number of fermionic Matsubara frequencies for G^2 calculations +g2_n_inu = 10 +# Number of Legendre coefficients for G^2 calculations +g2_n_l = 10 +# Block index combinations for G^2 calculations +g2_blocks = set([("up", "up"), ("up", "dn"), ("dn", "up")]) + +gf_struct = {"up" : orb_names, "dn" : orb_names} +print "Block structure of single-particle Green's functions:", gf_struct + +# Conversion from TRIQS to Pomerol notation for operator indices +# TRIQS: block_name, inner_index +# Pomerol: site_label, orbital_index, spin_name +index_converter = {} + +# Local degrees of freedom +index_converter.update({(sn, o) : ("loc", o, "down" if sn == "dn" else "up") + for sn, o in product(spin_names, orb_names)}) +# Bath degrees of freedom +index_converter.update({("B_" + sn, o) : ("bath", o, "down" if sn == "dn" else "up") + for sn, o in product(spin_names, orb_names)}) + +# Make PomerolED solver object +ed = PomerolED(index_converter, verbose = True) + +# Number of particles on the impurity +N = sum(n(sn, o) for sn, o in product(spin_names, orb_names)) + +# Local Hamiltonian +H_loc = h_int_kanamori(spin_names, orb_names, + np.array([[0, U-3*J], [U-3*J, 0]]), + np.array([[U, U-2*J], [U-2*J, U]]), + J, True) +H_loc -= mu*N + +# Bath Hamiltonian +H_bath = sum(epsilon[o] * n("B_" + sn, o) for sn, o in product(spin_names, orb_names)) + +# Hybridization Hamiltonian +H_hyb = sum( V[o1,o2] * c_dag("B_" + sn, o1) * c(sn, o2) + + np.conj(V[o2,o1]) * c_dag(sn, o1) * c("B_" + sn, o2) + for sn, o1, o2 in product(spin_names, orb_names, orb_names)) + +# Complete Hamiltonian +H = H_loc + H_hyb + H_bath + +# Diagonalize H +ed.diagonalize(H) + +# Compute G(i\omega) +G_iw = ed.G_iw(gf_struct, beta, n_iw) + +# Compute G(\tau) +G_tau = ed.G_tau(gf_struct, beta, n_tau) + +# Compute G(\omega) +G_w = ed.G_w(gf_struct, beta, energy_window, n_w, 0.01) + +########### +# G^{(2)} # +########### + +common_g2_params = {'gf_struct' : gf_struct, + 'beta' : beta, + 'blocks' : g2_blocks, + 'n_iw' : g2_n_iw} + +############################### +# G^{(2)}(i\omega;i\nu,i\nu') # +############################### + +# Compute G^{(2),ph}(i\omega;i\nu,i\nu'), AABB block order +G2_iw_inu_inup_ph_AABB = ed.G2_iw_inu_inup(channel = "PH", + block_order = "AABB", + n_inu = g2_n_inu, + **common_g2_params) + +# Compute G^{(2),ph}(i\omega;i\nu,i\nu'), ABBA block order +G2_iw_inu_inup_ph_ABBA = ed.G2_iw_inu_inup(channel = "PH", + block_order = "ABBA", + n_inu = g2_n_inu, + **common_g2_params) + +# Compute G^{(2),pp}(i\omega;i\nu,i\nu'), AABB block order +G2_iw_inu_inup_pp_AABB = ed.G2_iw_inu_inup(channel = "PP", + block_order = "AABB", + n_inu = g2_n_inu, + **common_g2_params) + +# Compute G^{(2),pp}(i\omega;i\nu,i\nu'), ABBA block order +G2_iw_inu_inup_pp_ABBA = ed.G2_iw_inu_inup(channel = "PP", + block_order = "ABBA", + n_inu = g2_n_inu, + **common_g2_params) + +######################### +# G^{(2)}(i\omega;l,l') # +######################### + +# Compute G^{(2),ph}(i\omega;l,l'), AABB block order +G2_iw_l_lp_ph_AABB = ed.G2_iw_l_lp(channel = "PH", + block_order = "AABB", + n_l = g2_n_l, + **common_g2_params) + +# Compute G^{(2),ph}(i\omega;l,l'), ABBA block order +G2_iw_l_lp_ph_ABBA = ed.G2_iw_l_lp(channel = "PH", + block_order = "ABBA", + n_l = g2_n_l, + **common_g2_params) + +# Compute G^{(2),pp}(i\omega;l,l'), AABB block order +G2_iw_l_lp_pp_AABB = ed.G2_iw_l_lp(channel = "PP", + block_order = "AABB", + n_l = g2_n_l, + **common_g2_params) + +# Compute G^{(2),pp}(i\omega;l,l'), ABBA block order +G2_iw_l_lp_pp_ABBA = ed.G2_iw_l_lp(channel = "PP", + block_order = "ABBA", + n_l = g2_n_l, + **common_g2_params) + +################ +# Save results # +################ + +if mpi.is_master_node(): + with HDFArchive('2band.h5', 'w') as ar: + ar['G_iw'] = G_iw + ar['G_tau'] = G_tau + ar['G_w'] = G_w + ar['G2_iw_inu_inup_ph_AABB'] = G2_iw_inu_inup_ph_AABB + ar['G2_iw_inu_inup_ph_ABBA'] = G2_iw_inu_inup_ph_ABBA + ar['G2_iw_inu_inup_pp_AABB'] = G2_iw_inu_inup_pp_AABB + ar['G2_iw_inu_inup_pp_ABBA'] = G2_iw_inu_inup_pp_ABBA + ar['G2_iw_l_lp_ph_AABB'] = G2_iw_l_lp_ph_AABB + ar['G2_iw_l_lp_ph_ABBA'] = G2_iw_l_lp_ph_ABBA + ar['G2_iw_l_lp_pp_AABB'] = G2_iw_l_lp_pp_AABB + ar['G2_iw_l_lp_pp_ABBA'] = G2_iw_l_lp_pp_ABBA diff --git a/example/anderson.py b/example/anderson.py new file mode 100644 index 0000000..31e3da8 --- /dev/null +++ b/example/anderson.py @@ -0,0 +1,140 @@ +from pytriqs.archive import HDFArchive +from pytriqs.gf import * +from pytriqs.operators import Operator, c, c_dag, n +from pytriqs.utility import mpi +from pytriqs.applications.impurity_solvers.pomerol2triqs import PomerolED +import numpy as np +from itertools import product + +# Single orbital Anderson model + +#################### +# Input parameters # +#################### + +beta = 10.0 # Inverse temperature +U = 2.0 # Coulomb repulsion +mu = 1.0 # Chemical potential + +# Levels of the bath sites +epsilon = [-1.0, 0, 1.0] +# Hopping amplitudes +V = [0.5, 0.5, 0.5] + +spin_names = ("up", "dn") + +# Number of Matsubara frequencies for GF calculation +n_iw = 1024 + +# Number of imaginary time slices for GF calculation +n_tau = 10001 + +# Energy window for real frequency GF calculation +energy_window = (-5, 5) +# Number of frequency points for real frequency GF calculation +n_w = 1000 + +# Number of bosonic Matsubara frequencies for G^2 calculations +g2_n_iw = 5 +# Number of fermionic Matsubara frequencies for G^2 calculations +g2_n_inu = 10 +# Number of Legendre coefficients for G^2 calculations +g2_n_l = 10 +# Block index combinations for G^2 calculations +g2_blocks = set([("up", "up"), ("up", "dn"), ("dn", "up")]) + +gf_struct = {'up' : [0], 'dn' : [0]} + +# Conversion from TRIQS to Pomerol notation for operator indices +# TRIQS: block_name, inner_index +# Pomerol: site_label, orbital_index, spin_name +index_converter = {} + +# Local degrees of freedom +index_converter.update({(sn, 0) : ("loc", 0, "down" if sn == "dn" else "up") for sn in spin_names}) +# Bath degrees of freedom +index_converter.update({("B%i_%s" % (k, sn), 0) : ("bath" + str(k), 0, "down" if sn == "dn" else "up") + for k, sn in product(range(len(epsilon)), spin_names)}) + +# Make PomerolED solver object +ed = PomerolED(index_converter, verbose = True) + +# Number of particles on the impurity +H_loc = -mu*(n('up', 0) + n('dn', 0)) + U * n('up', 0) * n('dn', 0) + +# Bath Hamiltonian +H_bath = sum(eps*n("B%i_%s" % (k, sn), 0) + for sn, (k, eps) in product(spin_names, enumerate(epsilon))) + +# Hybridization Hamiltonian +H_hyb = Operator() +for k, v in enumerate(V): + H_hyb += sum( v * c_dag("B%i_%s" % (k, sn), 0) * c(sn, 0) + + np.conj(v) * c_dag(sn, 0) * c("B%i_%s" % (k, sn), 0) + for sn in spin_names) + +# Complete Hamiltonian +H = H_loc + H_hyb + H_bath + +# Diagonalize H +ed.diagonalize(H) + +# Compute G(i\omega) +G_iw = ed.G_iw(gf_struct, beta, n_iw) + +# Compute G(\tau) +G_tau = ed.G_tau(gf_struct, beta, n_tau) + +# Compute G(\omega) +G_w = ed.G_w(gf_struct, beta, energy_window, n_w, 0.01) + +########### +# G^{(2)} # +########### + +common_g2_params = {'gf_struct' : gf_struct, + 'beta' : beta, + 'blocks' : g2_blocks, + 'n_iw' : g2_n_iw} + +############################### +# G^{(2)}(i\omega;i\nu,i\nu') # +############################### + +# Compute G^{(2),ph}(i\omega;i\nu,i\nu'), AABB block order +G2_iw_inu_inup_ph_AABB = ed.G2_iw_inu_inup(channel = "PH", + block_order = "AABB", + n_inu = g2_n_inu, + **common_g2_params) + +# Compute G^{(2),ph}(i\omega;i\nu,i\nu'), ABBA block order +G2_iw_inu_inup_ph_ABBA = ed.G2_iw_inu_inup(channel = "PH", + block_order = "ABBA", + n_inu = g2_n_inu, + **common_g2_params) + +# Compute G^{(2),pp}(i\omega;i\nu,i\nu'), AABB block order +G2_iw_inu_inup_pp_AABB = ed.G2_iw_inu_inup(channel = "PP", + block_order = "AABB", + n_inu = g2_n_inu, + **common_g2_params) + +# Compute G^{(2),pp}(i\omega;i\nu,i\nu'), ABBA block order +G2_iw_inu_inup_pp_ABBA = ed.G2_iw_inu_inup(channel = "PP", + block_order = "ABBA", + n_inu = g2_n_inu, + **common_g2_params) + +################ +# Save results # +################ + +if mpi.is_master_node(): + with HDFArchive('anderson.h5', 'w') as ar: + ar['G_iw'] = G_iw + ar['G_tau'] = G_tau + ar['G_w'] = G_w + ar['G2_iw_inu_inup_ph_AABB'] = G2_iw_inu_inup_ph_AABB + ar['G2_iw_inu_inup_ph_ABBA'] = G2_iw_inu_inup_ph_ABBA + ar['G2_iw_inu_inup_pp_AABB'] = G2_iw_inu_inup_pp_AABB + ar['G2_iw_inu_inup_pp_ABBA'] = G2_iw_inu_inup_pp_ABBA diff --git a/python/CMakeLists.txt b/python/CMakeLists.txt new file mode 100644 index 0000000..802f6a1 --- /dev/null +++ b/python/CMakeLists.txt @@ -0,0 +1,16 @@ +# where will the python end up in triqs? +set(python_destination pytriqs/applications/impurity_solvers/pomerol2triqs) + +# make the build_xxx, install python files... +triqs_prepare_local_pytriqs(${python_destination}) + +# sitecustomize.py for build +set(package_name "pytriqs.applications.impurity_solvers") +configure_file(${CMAKE_SOURCE_DIR}/cmake/sitecustomize.py ${CMAKE_CURRENT_BINARY_DIR}/sitecustomize.py @ONLY) + +# Build the python module for the solver +triqs_python_extension(pomerol2triqs ${python_destination}) +target_link_libraries(pomerol2triqs pomerol2triqs_c) +include_directories(${TRIQS_INCLUDE_ALL} ${CMAKE_CURRENT_SOURCE_DIR} ${pomerol_INCLUDE_DIRS}) +triqs_set_rpath_for_target(pomerol2triqs) +install(TARGETS pomerol2triqs DESTINATION ${TRIQS_PYTHON_LIB_DEST_ROOT}/${python_destination}) diff --git a/python/__init__.py b/python/__init__.py new file mode 100644 index 0000000..f4bc587 --- /dev/null +++ b/python/__init__.py @@ -0,0 +1,7 @@ +r""" +DOC +""" + +from pomerol2triqs import PomerolED + +__all__ = ['PomerolED'] diff --git a/python/pomerol2triqs_converters.hxx b/python/pomerol2triqs_converters.hxx new file mode 100644 index 0000000..17edfce --- /dev/null +++ b/python/pomerol2triqs_converters.hxx @@ -0,0 +1,219 @@ +// DO NOT EDIT +// Generated automatically using libclang using the command : +// c++2py.py ../c++/pomerol_ed.hpp -I../../../pomerol/installed/include -I/usr/include/eigen3 -I../c++ --only_converters -p -mpytriqs.applications.impurity_solvers.pomerol2triqs -o pomerol2triqs --moduledoc "TRIQS wrapper around Pomerol ED library" + + +// --- C++ Python converter for g2_iw_l_lp_params_t +#include +#include +#include + +namespace triqs { namespace py_tools { + +template <> struct py_converter { + static PyObject *c2py(g2_iw_l_lp_params_t const & x) { + PyObject * d = PyDict_New(); + PyDict_SetItemString( d, "gf_struct" , convert_to_python(x.gf_struct)); + PyDict_SetItemString( d, "beta" , convert_to_python(x.beta)); + PyDict_SetItemString( d, "channel" , convert_to_python(x.channel)); + PyDict_SetItemString( d, "block_order", convert_to_python(x.block_order)); + PyDict_SetItemString( d, "blocks" , convert_to_python(x.blocks)); + PyDict_SetItemString( d, "n_iw" , convert_to_python(x.n_iw)); + PyDict_SetItemString( d, "n_l" , convert_to_python(x.n_l)); + PyDict_SetItemString( d, "n_inu_sum" , convert_to_python(x.n_inu_sum)); + PyDict_SetItemString( d, "inu_sum_tol", convert_to_python(x.inu_sum_tol)); + return d; + } + + template static void _get_optional(PyObject *dic, const char *name, T &r, U const &init_default) { + if (PyDict_Contains(dic, pyref::string(name))) + r = convert_from_python(PyDict_GetItemString(dic, name)); + else + r = init_default; + } + + template static void _get_optional(PyObject *dic, const char *name, T &r) { + if (PyDict_Contains(dic, pyref::string(name))) + r = convert_from_python(PyDict_GetItemString(dic, name)); + else + r = T{}; + } + + static g2_iw_l_lp_params_t py2c(PyObject *dic) { + g2_iw_l_lp_params_t res; + res.gf_struct = convert_from_python(PyDict_GetItemString(dic, "gf_struct")); + res.beta = convert_from_python(PyDict_GetItemString(dic, "beta")); + _get_optional(dic, "channel" , res.channel ,PH); + _get_optional(dic, "block_order", res.block_order ,AABB); + _get_optional(dic, "blocks" , res.blocks ,g2_blocks_t{}); + _get_optional(dic, "n_iw" , res.n_iw ,30); + _get_optional(dic, "n_l" , res.n_l ,20); + _get_optional(dic, "n_inu_sum" , res.n_inu_sum ,500); + _get_optional(dic, "inu_sum_tol", res.inu_sum_tol ,1e-6); + return res; + } + + template + static void _check(PyObject *dic, std::stringstream &fs, int &err, const char *name, const char *tname) { + if (!convertible_from_python(PyDict_GetItemString(dic, name), false)) + fs << "\n" << ++err << " The parameter " << name << " does not have the right type : expecting " << tname + << " in C++, but got '" << PyDict_GetItemString(dic, name)->ob_type->tp_name << "' in Python."; + } + + template + static void _check_mandatory(PyObject *dic, std::stringstream &fs, int &err, const char *name, const char *tname) { + if (!PyDict_Contains(dic, pyref::string(name))) + fs << "\n" << ++err << " Mandatory parameter " << name << " is missing."; + else _check(dic,fs,err,name,tname); + } + + template + static void _check_optional(PyObject *dic, std::stringstream &fs, int &err, const char *name, const char *tname) { + if (PyDict_Contains(dic, pyref::string(name))) _check(dic, fs, err, name, tname); + } + + static bool is_convertible(PyObject *dic, bool raise_exception) { + if (dic == nullptr or !PyDict_Check(dic)) { + if (raise_exception) { PyErr_SetString(PyExc_TypeError, "The function must be called with named arguments");} + return false; + } + std::stringstream fs, fs2; int err=0; + +#ifndef TRIQS_ALLOW_UNUSED_PARAMETERS + std::vector ks, all_keys = {"gf_struct","beta","channel","block_order","blocks","n_iw","n_l","n_inu_sum","inu_sum_tol"}; + pyref keys = PyDict_Keys(dic); + if (!convertible_from_python>(keys, true)) { + fs << "\nThe dict keys are not strings"; + goto _error; + } + ks = convert_from_python>(keys); + for (auto & k : ks) + if (std::find(all_keys.begin(), all_keys.end(), k) == all_keys.end()) + fs << "\n"<< ++err << " The parameter '" << k << "' is not recognized."; +#endif + + _check_mandatory(dic, fs, err, "gf_struct" , "gf_struct_t"); + _check_mandatory(dic, fs, err, "beta" , "double"); + _check_optional (dic, fs, err, "channel" , "pomerol2triqs::channel_t"); + _check_optional (dic, fs, err, "block_order", "pomerol2triqs::block_order_t"); + _check_optional (dic, fs, err, "blocks" , "g2_blocks_t"); + _check_optional (dic, fs, err, "n_iw" , "int"); + _check_optional (dic, fs, err, "n_l" , "int"); + _check_optional (dic, fs, err, "n_inu_sum" , "int"); + _check_optional (dic, fs, err, "inu_sum_tol", "double"); + if (err) goto _error; + return true; + + _error: + fs2 << "\n---- There " << (err > 1 ? "are " : "is ") << err<< " error"<<(err >1 ?"s" : "")<< " in Python -> C++ transcription for the class g2_iw_l_lp_params_t\n" < +#include +#include + +namespace triqs { namespace py_tools { + +template <> struct py_converter { + static PyObject *c2py(g2_iw_inu_inup_params_t const & x) { + PyObject * d = PyDict_New(); + PyDict_SetItemString( d, "gf_struct" , convert_to_python(x.gf_struct)); + PyDict_SetItemString( d, "beta" , convert_to_python(x.beta)); + PyDict_SetItemString( d, "channel" , convert_to_python(x.channel)); + PyDict_SetItemString( d, "block_order", convert_to_python(x.block_order)); + PyDict_SetItemString( d, "blocks" , convert_to_python(x.blocks)); + PyDict_SetItemString( d, "n_iw" , convert_to_python(x.n_iw)); + PyDict_SetItemString( d, "n_inu" , convert_to_python(x.n_inu)); + return d; + } + + template static void _get_optional(PyObject *dic, const char *name, T &r, U const &init_default) { + if (PyDict_Contains(dic, pyref::string(name))) + r = convert_from_python(PyDict_GetItemString(dic, name)); + else + r = init_default; + } + + template static void _get_optional(PyObject *dic, const char *name, T &r) { + if (PyDict_Contains(dic, pyref::string(name))) + r = convert_from_python(PyDict_GetItemString(dic, name)); + else + r = T{}; + } + + static g2_iw_inu_inup_params_t py2c(PyObject *dic) { + g2_iw_inu_inup_params_t res; + res.gf_struct = convert_from_python(PyDict_GetItemString(dic, "gf_struct")); + res.beta = convert_from_python(PyDict_GetItemString(dic, "beta")); + _get_optional(dic, "channel" , res.channel ,PH); + _get_optional(dic, "block_order", res.block_order ,AABB); + _get_optional(dic, "blocks" , res.blocks ,g2_blocks_t{}); + _get_optional(dic, "n_iw" , res.n_iw ,30); + _get_optional(dic, "n_inu" , res.n_inu ,30); + return res; + } + + template + static void _check(PyObject *dic, std::stringstream &fs, int &err, const char *name, const char *tname) { + if (!convertible_from_python(PyDict_GetItemString(dic, name), false)) + fs << "\n" << ++err << " The parameter " << name << " does not have the right type : expecting " << tname + << " in C++, but got '" << PyDict_GetItemString(dic, name)->ob_type->tp_name << "' in Python."; + } + + template + static void _check_mandatory(PyObject *dic, std::stringstream &fs, int &err, const char *name, const char *tname) { + if (!PyDict_Contains(dic, pyref::string(name))) + fs << "\n" << ++err << " Mandatory parameter " << name << " is missing."; + else _check(dic,fs,err,name,tname); + } + + template + static void _check_optional(PyObject *dic, std::stringstream &fs, int &err, const char *name, const char *tname) { + if (PyDict_Contains(dic, pyref::string(name))) _check(dic, fs, err, name, tname); + } + + static bool is_convertible(PyObject *dic, bool raise_exception) { + if (dic == nullptr or !PyDict_Check(dic)) { + if (raise_exception) { PyErr_SetString(PyExc_TypeError, "The function must be called with named arguments");} + return false; + } + std::stringstream fs, fs2; int err=0; + +#ifndef TRIQS_ALLOW_UNUSED_PARAMETERS + std::vector ks, all_keys = {"gf_struct","beta","channel","block_order","blocks","n_iw","n_inu"}; + pyref keys = PyDict_Keys(dic); + if (!convertible_from_python>(keys, true)) { + fs << "\nThe dict keys are not strings"; + goto _error; + } + ks = convert_from_python>(keys); + for (auto & k : ks) + if (std::find(all_keys.begin(), all_keys.end(), k) == all_keys.end()) + fs << "\n"<< ++err << " The parameter '" << k << "' is not recognized."; +#endif + + _check_mandatory(dic, fs, err, "gf_struct" , "gf_struct_t"); + _check_mandatory(dic, fs, err, "beta" , "double"); + _check_optional (dic, fs, err, "channel" , "pomerol2triqs::channel_t"); + _check_optional (dic, fs, err, "block_order", "pomerol2triqs::block_order_t"); + _check_optional (dic, fs, err, "blocks" , "g2_blocks_t"); + _check_optional (dic, fs, err, "n_iw" , "int"); + _check_optional (dic, fs, err, "n_inu" , "int"); + if (err) goto _error; + return true; + + _error: + fs2 << "\n---- There " << (err > 1 ? "are " : "is ") << err<< " error"<<(err >1 ?"s" : "")<< " in Python -> C++ transcription for the class g2_iw_inu_inup_params_t\n" < +#include +#include +#include +#include +#include +#include +#include +#include +#include "./pomerol2triqs_converters.hxx" +""") +module.add_using("namespace pomerol2triqs") +module.add_using("namespace Pomerol") + +module.add_enum("spin", ["down", "up"], "Pomerol", "Spin projection") +module.add_enum("block_order_t", ["AABB", "ABBA"], "pomerol2triqs", "G^{(2)} block order") +module.add_enum("channel_t", ["PP", "PH", "AllFermionic"], "pomerol2triqs", "G^{(2)} channel") + +# The class pomerol_ed +c = class_( + py_type = "PomerolED", # name of the python class + c_type = "pomerol_ed", # name of the C++ class + doc = r"Main solver class of pomerol2triqs", # doc of the C++ class +) + +c.add_constructor("""(index_converter_t index_converter, bool verbose = false)""", + doc = r"""Create a new solver object""") + +c.add_method("""void diagonalize(many_body_op_t hamiltonian, bool ignore_symmetries = false)""", + doc = r"""Diagonalize Hamiltonian""") + +c.add_method("""void diagonalize(many_body_op_t hamiltonian, std::vector integrals_of_motion)""", + doc = r"""Diagonalize Hamiltonian""") + +c.add_method("""block_gf G_iw (gf_struct_t gf_struct, double beta, int n_iw)""", + doc = r"""Green's function in Matsubara frequencies""") + +c.add_method("""block_gf G_tau (gf_struct_t gf_struct, double beta, int n_tau)""", + doc = r"""Green's function in imaginary time""") + +c.add_method("""block_gf G_w(gf_struct_t gf_struct, double beta, std::pair energy_window, int n_w, double im_shift = 0)""", + doc = r"""Retarded Green's function on real energy axis""") + +c.add_method("""block2_gf, tensor_valued<4>> G2_iw_inu_inup(**pomerol2triqs::g2_iw_inu_inup_params_t)""", + doc = r"""Two-particle Green's function, Matsubara frequencies""") + +c.add_method("""block2_gf, tensor_valued<4>> G2_iw_l_lp(**pomerol2triqs::g2_iw_l_lp_params_t)""", + doc = r"""Two-particle Green's function, bosonic Matsubara frequency + Legendre coefficients""") + +module.add_class(c) + +module.generate_code() diff --git a/python/pomerol2triqs_parameters.rst b/python/pomerol2triqs_parameters.rst new file mode 100644 index 0000000..d77e38c --- /dev/null +++ b/python/pomerol2triqs_parameters.rst @@ -0,0 +1,41 @@ ++----------------+------------------------------+--------------------+------------------------------------------------------------------+ +| Parameter Name | Type | Default | Documentation | ++================+==============================+====================+==================================================================+ +| gf_struct | gf_struct_t | -- | Structure of G^2 blocks. | ++----------------+------------------------------+--------------------+------------------------------------------------------------------+ +| beta | double | -- | Inverse temperature | ++----------------+------------------------------+--------------------+------------------------------------------------------------------+ +| channel | channel_t | PH | Channel in which Matsubara frequency representation is defined. | ++----------------+------------------------------+--------------------+------------------------------------------------------------------+ +| block_order | block_order_t | AABB | Order of block indices in the definition of G^2. | ++----------------+------------------------------+--------------------+------------------------------------------------------------------+ +| blocks | g2_blocks_t | measure all blocks | List of block index pairs of G^2 to measure. | ++----------------+------------------------------+--------------------+------------------------------------------------------------------+ +| n_iw | int | 30 | Number of bosonic Matsubara frequencies. | ++----------------+------------------------------+--------------------+------------------------------------------------------------------+ +| n_l | int | 20 | Number of Legendre coefficients. | ++----------------+------------------------------+--------------------+------------------------------------------------------------------+ +| n_inu_sum | int | 500 | Maximum number of positive Matsubara frequencies in summation. | ++----------------+------------------------------+--------------------+------------------------------------------------------------------+ +| inu_sum_tol | double | 1e-6 | Tolerance for Matsubara frequency summation. | ++----------------+------------------------------+--------------------+------------------------------------------------------------------+ + + + ++----------------+------------------------------+--------------------+------------------------------------------------------------------+ +| Parameter Name | Type | Default | Documentation | ++================+==============================+====================+==================================================================+ +| gf_struct | gf_struct_t | -- | Structure of G^2 blocks. | ++----------------+------------------------------+--------------------+------------------------------------------------------------------+ +| beta | double | -- | Inverse temperature | ++----------------+------------------------------+--------------------+------------------------------------------------------------------+ +| channel | channel_t | PH | Channel in which Matsubara frequency representation is defined. | ++----------------+------------------------------+--------------------+------------------------------------------------------------------+ +| block_order | block_order_t | AABB | Order of block indices in the definition of G^2. | ++----------------+------------------------------+--------------------+------------------------------------------------------------------+ +| blocks | g2_blocks_t | measure all blocks | List of block index pairs of G^2 to measure. | ++----------------+------------------------------+--------------------+------------------------------------------------------------------+ +| n_iw | int | 30 | Number of bosonic Matsubara frequencies. | ++----------------+------------------------------+--------------------+------------------------------------------------------------------+ +| n_inu | int | 30 | Number of fermionic Matsubara frequencies. | ++----------------+------------------------------+--------------------+------------------------------------------------------------------+ diff --git a/test/CMakeLists.txt b/test/CMakeLists.txt new file mode 100644 index 0000000..8e5f91a --- /dev/null +++ b/test/CMakeLists.txt @@ -0,0 +1 @@ +add_subdirectory(python) diff --git a/test/python/CMakeLists.txt b/test/python/CMakeLists.txt new file mode 100644 index 0000000..de8460f --- /dev/null +++ b/test/python/CMakeLists.txt @@ -0,0 +1,13 @@ +find_package(TriqsTest) + +# Copy h5 files to binary dir +FILE(GLOB all_h5_files RELATIVE ${CMAKE_CURRENT_SOURCE_DIR} *.h5) +file(COPY ${CMAKE_CURRENT_SOURCE_DIR}/${all_h5_files} DESTINATION ${CMAKE_CURRENT_BINARY_DIR}) + +triqs_add_python_test(anderson_gf) +triqs_add_python_test(slater_gf) +triqs_add_python_test(wick) + +foreach(TEST_MPI_NUMPROC 1 2 4) + triqs_add_python_test(anderson_g2_matsubara) +endforeach() diff --git a/test/python/anderson_g2_matsubara.py b/test/python/anderson_g2_matsubara.py new file mode 100644 index 0000000..3a0814f --- /dev/null +++ b/test/python/anderson_g2_matsubara.py @@ -0,0 +1,114 @@ +from pytriqs.archive import HDFArchive +from pytriqs.gf import * +from pytriqs.operators import Operator, c, c_dag, n +from pytriqs.utility import mpi +from pytriqs.applications.impurity_solvers.pomerol2triqs import PomerolED +from pytriqs.utility.comparison_tests import * +import numpy as np +from itertools import product + +# Single orbital Anderson model + +#################### +# Input parameters # +#################### + +beta = 10.0 # Inverse temperature +U = 2.0 # Coulomb repulsion +mu = 1.0 # Chemical potential + +# Levels of the bath sites +epsilon = [-1.0, 1.0] +# Hopping amplitudes +V = [0.5, 0.5] + +spin_names = ("up", "dn") + +# Number of bosonic Matsubara frequencies for G^2 calculations +g2_n_iw = 5 +# Number of fermionic Matsubara frequencies for G^2 calculations +g2_n_inu = 5 +# Block index combinations for G^2 calculations +g2_blocks = set([("up", "up"), ("dn", "dn"), ("up", "dn")]) + +# GF structure +gf_struct = {'up' : [0], 'dn' : [0]} + +# Conversion from TRIQS to Pomerol notation for operator indices +index_converter = {} +index_converter.update({(sn, 0) : ("loc", 0, "down" if sn == "dn" else "up") for sn in spin_names}) +index_converter.update({("B%i_%s" % (k, sn), 0) : ("bath" + str(k), 0, "down" if sn == "dn" else "up") + for k, sn in product(range(len(epsilon)), spin_names)}) + +# Make PomerolED solver object +ed = PomerolED(index_converter, verbose = True) + +# Number of particles on the impurity +H_loc = -mu*(n('up', 0) + n('dn', 0)) + U * n('up', 0) * n('dn', 0) + +# Bath Hamiltonian +H_bath = sum(eps*n("B%i_%s" % (k, sn), 0) + for sn, (k, eps) in product(spin_names, enumerate(epsilon))) + +# Hybridization Hamiltonian +H_hyb = Operator() +for k, v in enumerate(V): + H_hyb += sum( v * c_dag("B%i_%s" % (k, sn), 0) * c(sn, 0) + + np.conj(v) * c_dag(sn, 0) * c("B%i_%s" % (k, sn), 0) + for sn in spin_names) + +# Complete Hamiltonian +H = H_loc + H_hyb + H_bath + +# Diagonalize H +ed.diagonalize(H) + +########### +# G^{(2)} # +########### + +common_g2_params = {'gf_struct' : gf_struct, + 'beta' : beta, + 'blocks' : g2_blocks, + 'n_iw' : g2_n_iw} + +# Compute G^{(2),ph}(i\omega;i\nu,i\nu'), AABB block order +G2_iw_inu_inup_ph_AABB = ed.G2_iw_inu_inup(channel = "PH", + block_order = "AABB", + n_inu = g2_n_inu, + **common_g2_params) + +# Compute G^{(2),ph}(i\omega;i\nu,i\nu'), ABBA block order +G2_iw_inu_inup_ph_ABBA = ed.G2_iw_inu_inup(channel = "PH", + block_order = "ABBA", + n_inu = g2_n_inu, + **common_g2_params) + +# Compute G^{(2),pp}(i\omega;i\nu,i\nu'), AABB block order +G2_iw_inu_inup_pp_AABB = ed.G2_iw_inu_inup(channel = "PP", + block_order = "AABB", + n_inu = g2_n_inu, + **common_g2_params) + +# Compute G^{(2),pp}(i\omega;i\nu,i\nu'), ABBA block order +G2_iw_inu_inup_pp_ABBA = ed.G2_iw_inu_inup(channel = "PP", + block_order = "ABBA", + n_inu = g2_n_inu, + **common_g2_params) + +if mpi.is_master_node(): + with HDFArchive('anderson_g2_matsubara_np%i.out.h5' % mpi.size, 'w') as ar: + ar['H'] = H + ar['G2_iw_inu_inup_ph_AABB'] = G2_iw_inu_inup_ph_AABB + ar['G2_iw_inu_inup_ph_ABBA'] = G2_iw_inu_inup_ph_ABBA + ar['G2_iw_inu_inup_pp_AABB'] = G2_iw_inu_inup_pp_AABB + ar['G2_iw_inu_inup_pp_ABBA'] = G2_iw_inu_inup_pp_ABBA + + with HDFArchive("anderson_g2_matsubara.ref.h5", 'r') as ar: + assert (ar['H'] - H).is_zero() + for bn1, bn2 in g2_blocks: + assert_gfs_are_close(ar['G2_iw_inu_inup_ph_AABB'][bn1, bn2], G2_iw_inu_inup_ph_AABB[bn1, bn2]) + assert_gfs_are_close(ar['G2_iw_inu_inup_ph_ABBA'][bn1, bn2], G2_iw_inu_inup_ph_ABBA[bn1, bn2]) + assert_gfs_are_close(ar['G2_iw_inu_inup_pp_AABB'][bn1, bn2], G2_iw_inu_inup_pp_AABB[bn1, bn2]) + assert_gfs_are_close(ar['G2_iw_inu_inup_pp_ABBA'][bn1, bn2], G2_iw_inu_inup_pp_ABBA[bn1, bn2]) + diff --git a/test/python/anderson_g2_matsubara.ref.h5 b/test/python/anderson_g2_matsubara.ref.h5 new file mode 100644 index 0000000..df2845f Binary files /dev/null and b/test/python/anderson_g2_matsubara.ref.h5 differ diff --git a/test/python/anderson_gf.py b/test/python/anderson_gf.py new file mode 100644 index 0000000..fd6cad6 --- /dev/null +++ b/test/python/anderson_gf.py @@ -0,0 +1,90 @@ +from pytriqs.archive import HDFArchive +from pytriqs.gf import * +from pytriqs.operators import Operator, c, c_dag, n +from pytriqs.utility import mpi +from pytriqs.applications.impurity_solvers.pomerol2triqs import PomerolED +from pytriqs.utility.comparison_tests import * +import numpy as np +from itertools import product + +# Single orbital Anderson model + +#################### +# Input parameters # +#################### + +beta = 10.0 # Inverse temperature +U = 2.0 # Coulomb repulsion +mu = 1.0 # Chemical potential + +# Levels of the bath sites +epsilon = [-1.0, 1.0] +# Hopping amplitudes +V = [0.5, 0.5] + +spin_names = ("up", "dn") + +# Number of Matsubara frequencies for GF calculation +n_iw = 200 + +# Number of imaginary time slices for GF calculation +n_tau = 1001 + +# Energy window for real frequency GF calculation +energy_window = (-5, 5) +# Number of frequency points for real frequency GF calculation +n_w = 1000 + +# GF structure +gf_struct = {'up' : [0], 'dn' : [0]} + +# Conversion from TRIQS to Pomerol notation for operator indices +index_converter = {} +index_converter.update({(sn, 0) : ("loc", 0, "down" if sn == "dn" else "up") for sn in spin_names}) +index_converter.update({("B%i_%s" % (k, sn), 0) : ("bath" + str(k), 0, "down" if sn == "dn" else "up") + for k, sn in product(range(len(epsilon)), spin_names)}) + +# Make PomerolED solver object +ed = PomerolED(index_converter, verbose = True) + +# Number of particles on the impurity +H_loc = -mu*(n('up', 0) + n('dn', 0)) + U * n('up', 0) * n('dn', 0) + +# Bath Hamiltonian +H_bath = sum(eps*n("B%i_%s" % (k, sn), 0) + for sn, (k, eps) in product(spin_names, enumerate(epsilon))) + +# Hybridization Hamiltonian +H_hyb = Operator() +for k, v in enumerate(V): + H_hyb += sum( v * c_dag("B%i_%s" % (k, sn), 0) * c(sn, 0) + + np.conj(v) * c_dag(sn, 0) * c("B%i_%s" % (k, sn), 0) + for sn in spin_names) + +# Complete Hamiltonian +H = H_loc + H_hyb + H_bath + +# Diagonalize H +ed.diagonalize(H) + +# Compute G(i\omega) +G_iw = ed.G_iw(gf_struct, beta, n_iw) + +# Compute G(\tau) +G_tau = ed.G_tau(gf_struct, beta, n_tau) + +# Compute G(\omega) +G_w = ed.G_w(gf_struct, beta, energy_window, n_w, 0.01) + +if mpi.is_master_node(): + with HDFArchive('anderson_gf.out.h5', 'w') as ar: + ar['H'] = H + ar['G_iw'] = G_iw + ar['G_tau'] = G_tau + ar['G_w'] = G_w + + with HDFArchive("anderson_gf.ref.h5", 'r') as ar: + assert (ar['H'] - H).is_zero() + assert_block_gfs_are_close(ar['G_iw'], G_iw) + assert_block_gfs_are_close(ar['G_tau'], G_tau) + assert_block_gfs_are_close(ar['G_w'], G_w) diff --git a/test/python/anderson_gf.ref.h5 b/test/python/anderson_gf.ref.h5 new file mode 100644 index 0000000..f8cfbfe Binary files /dev/null and b/test/python/anderson_gf.ref.h5 differ diff --git a/test/python/slater_gf.py b/test/python/slater_gf.py new file mode 100644 index 0000000..3f3e186 --- /dev/null +++ b/test/python/slater_gf.py @@ -0,0 +1,96 @@ +from pytriqs.archive import HDFArchive +from pytriqs.gf import * +from pytriqs.operators import * +from pytriqs.operators.util.op_struct import set_operator_structure, get_mkind +from pytriqs.operators.util.U_matrix import U_matrix +from pytriqs.operators.util.hamiltonians import h_int_slater +from pytriqs.operators.util.observables import N_op, S_op, L_op +from pytriqs.applications.impurity_solvers.pomerol2triqs import PomerolED +from pytriqs.utility import mpi +from pytriqs.utility.comparison_tests import * +from itertools import product + +# 5-orbital atom with Slater interaction term + +#################### +# Input parameters # +#################### + +L = 2 # Angular momentum +beta = 10.0 # Inverse temperature +mu = 32.5 # Chemical potential (3 electrons in 5 bands) +# Slater parameters +U = 5.0 +J = 0.1 +F0 = U +F2 = J*(14.0/(1.0 + 0.625)) +F4 = F2*0.625 + +spin_names = ("up", "dn") +orb_names = range(-L, L+1) +U_mat = U_matrix(L, radial_integrals = [F0,F2,F4], basis='spherical') + +# Number of Matsubara frequencies for GF calculation +n_iw = 200 + +# Number of imaginary time slices for GF calculation +n_tau = 1001 + +# Energy window for real frequency GF calculation +energy_window = (-5, 5) +# Number of frequency points for real frequency GF calculation +n_w = 1000 + +# GF structure +gf_struct = set_operator_structure(spin_names, orb_names, False) + +mkind = get_mkind(False, None) +# Conversion from TRIQS to Pomerol notation for operator indices +index_converter = {mkind(sn, bn) : ("atom", bi, "down" if sn == "dn" else "up") + for sn, (bi, bn) in product(spin_names, enumerate(orb_names))} + +# Make PomerolED solver object +ed = PomerolED(index_converter, verbose = True) + +# Hamiltonian +H = h_int_slater(spin_names, orb_names, U_mat, False) + +# Number of particles +N = N_op(spin_names, orb_names, False) + +# z-component of spin +Sz = S_op('z', spin_names, orb_names, False) + +# z-component of angular momentum +Lz = L_op('z', spin_names, orb_names, off_diag = False, basis = 'spherical') + +# Double check that we are actually using integrals of motion +h_comm = lambda op: H*op - op*H +assert h_comm(N).is_zero() +assert h_comm(Sz).is_zero() +assert h_comm(Lz).is_zero() + +# Diagonalize H + +# Do not split H into blocks (uncomment to generate reference data) +#ed.diagonalize(H, True) + +ed.diagonalize(H, [N, Sz, Lz]) + +# Compute G(i\omega) +G_iw = ed.G_iw(gf_struct, beta, n_iw) + +# Compute G(\tau) +G_tau = ed.G_tau(gf_struct, beta, n_tau) + +if mpi.is_master_node(): + with HDFArchive('slater_gf.out.h5', 'w') as ar: + ar['H'] = H + ar['G_iw'] = G_iw + ar['G_tau'] = G_tau + + with HDFArchive("slater_gf.ref.h5", 'r') as ar: + assert (ar['H'] - H).is_zero() + assert_block_gfs_are_close(ar['G_iw'], G_iw) + assert_block_gfs_are_close(ar['G_tau'], G_tau) + diff --git a/test/python/slater_gf.ref.h5 b/test/python/slater_gf.ref.h5 new file mode 100644 index 0000000..5669aa4 Binary files /dev/null and b/test/python/slater_gf.ref.h5 differ diff --git a/test/python/wick.py b/test/python/wick.py new file mode 100644 index 0000000..05df5c1 --- /dev/null +++ b/test/python/wick.py @@ -0,0 +1,148 @@ +from pytriqs.archive import HDFArchive +from pytriqs.gf import * +from pytriqs.operators import Operator, c, c_dag, n +from pytriqs.utility import mpi +from pytriqs.applications.impurity_solvers.pomerol2triqs import PomerolED +from pytriqs.utility.comparison_tests import * +import numpy as np +from itertools import product + +# Non-interacting impurity in a bath + +#################### +# Input parameters # +#################### + +beta = 10.0 # Inverse temperature +e_d = 0.1 # Local level +h = 0.15 # Magnetic field + +# Levels of the bath sites +epsilon = [-1.0, 1.0] +# Hopping amplitudes +V = [0.5, 0.5] + +spin_names = ("up", "dn") + +# Number of Matsubara frequencies for GF calculation +n_iw = 200 + +# Number of bosonic Matsubara frequencies for G^2 calculations + +g2_n_iw = 5 + +# Number of fermionic Matsubara frequencies for G^2 calculations +g2_n_inu = 5 + +# GF structure +gf_struct = {'up' : [0], 'dn' : [0]} + +# Conversion from TRIQS to Pomerol notation for operator indices +index_converter = {} +index_converter.update({(sn, 0) : ("loc", 0, "down" if sn == "dn" else "up") for sn in spin_names}) +index_converter.update({("B%i_%s" % (k, sn), 0) : ("bath" + str(k), 0, "down" if sn == "dn" else "up") + for k, sn in product(range(len(epsilon)), spin_names)}) + +# Make PomerolED solver object +ed = PomerolED(index_converter, verbose = True) + +# Local Hamiltonian +H_loc = e_d*(n('up', 0) + n('dn', 0)) + h*(n('up', 0) - n('dn', 0)) + +# Bath Hamiltonian +H_bath = sum(eps*n("B%i_%s" % (k, sn), 0) + for sn, (k, eps) in product(spin_names, enumerate(epsilon))) + +# Hybridization Hamiltonian +H_hyb = Operator() +for k, v in enumerate(V): + H_hyb += sum( v * c_dag("B%i_%s" % (k, sn), 0) * c(sn, 0) + + np.conj(v) * c_dag(sn, 0) * c("B%i_%s" % (k, sn), 0) + for sn in spin_names) + +# Complete Hamiltonian +H = H_loc + H_hyb + H_bath + +# Diagonalize H +ed.diagonalize(H) + +# Compute G(i\omega) +G_iw = ed.G_iw(gf_struct, beta, n_iw) + +########### +# G^{(2)} # +########### + +common_g2_params = {'gf_struct' : gf_struct, + 'beta' : beta, + 'n_iw' : g2_n_iw} + +# Compute G^{(2),ph}(i\omega;i\nu,i\nu'), AABB block order +G2_iw_inu_inup_ph_AABB = ed.G2_iw_inu_inup(channel = "PH", + block_order = "AABB", + n_inu = g2_n_inu, + **common_g2_params) + +# Compute G^{(2),ph}(i\omega;i\nu,i\nu'), ABBA block order +G2_iw_inu_inup_ph_ABBA = ed.G2_iw_inu_inup(channel = "PH", + block_order = "ABBA", + n_inu = g2_n_inu, + **common_g2_params) + +# Compute G^{(2),pp}(i\omega;i\nu,i\nu'), AABB block order +G2_iw_inu_inup_pp_AABB = ed.G2_iw_inu_inup(channel = "PP", + block_order = "AABB", + n_inu = g2_n_inu, + **common_g2_params) + +# Compute G^{(2),pp}(i\omega;i\nu,i\nu'), ABBA block order +G2_iw_inu_inup_pp_ABBA = ed.G2_iw_inu_inup(channel = "PP", + block_order = "ABBA", + n_inu = g2_n_inu, + **common_g2_params) + +b_mesh = MeshImFreq(beta = beta, S = 'Boson', n_max = g2_n_iw) +f_mesh = MeshImFreq(beta = beta, S = 'Fermion', n_max = g2_n_inu) +g2_mesh = MeshProduct(b_mesh, f_mesh, f_mesh) + +# Check meshes of G2 containers +for bn1, bn2 in product(spin_names, spin_names): + assert G2_iw_inu_inup_ph_AABB[bn1, bn2].mesh == g2_mesh + assert G2_iw_inu_inup_ph_ABBA[bn1, bn2].mesh == g2_mesh + assert G2_iw_inu_inup_pp_AABB[bn1, bn2].mesh == g2_mesh + assert G2_iw_inu_inup_pp_ABBA[bn1, bn2].mesh == g2_mesh + +G2_iw_inu_inup_ph_AABB_wick = G2_iw_inu_inup_ph_AABB.copy() +G2_iw_inu_inup_ph_ABBA_wick = G2_iw_inu_inup_ph_AABB.copy() +G2_iw_inu_inup_pp_AABB_wick = G2_iw_inu_inup_ph_AABB.copy() +G2_iw_inu_inup_pp_ABBA_wick = G2_iw_inu_inup_ph_AABB.copy() + +G = lambda s, i: G_iw[s].data[i + n_iw, 0, 0] + +g2_mesh_enum = product(enumerate(b_mesh), enumerate(f_mesh), enumerate(f_mesh)) +for (m, w), (i, nu), (ip, nup) in g2_mesh_enum: + m -= g2_n_iw - 1 + i -= g2_n_inu + ip -= g2_n_inu + d_w_0 = int(m == 0) + d_w_nu_nup = int(m == i + ip + 1) + d_nu_nup = int(i == ip) + + + for s1, s2 in product(spin_names, spin_names): + d_s1_s2 = int(s1 == s2) + + G2_iw_inu_inup_ph_AABB_wick[s1, s2].data[m + (g2_n_iw - 1), i + g2_n_inu, ip + g2_n_inu] = \ + beta * d_w_0 * G(s1, i) * G(s2, ip) - beta * d_nu_nup * d_s1_s2 * G(s1, m + i) * G(s2, i) + G2_iw_inu_inup_ph_ABBA_wick[s1, s2].data[m + (g2_n_iw - 1), i + g2_n_inu, ip + g2_n_inu] = \ + beta * d_w_0 * d_s1_s2 * G(s1, i) * G(s2, ip) - beta * d_nu_nup * G(s2, m + i) * G(s1, i) + G2_iw_inu_inup_pp_AABB_wick[s1, s2].data[m + (g2_n_iw - 1), i + g2_n_inu, ip + g2_n_inu] = \ + beta * d_w_nu_nup * G(s1, i) * G(s2, ip) - beta * d_nu_nup * d_s1_s2 * G(s1, m - i - 1) * G(s2, i) + G2_iw_inu_inup_pp_ABBA_wick[s1, s2].data[m + (g2_n_iw - 1), i + g2_n_inu, ip + g2_n_inu] = \ + beta * d_w_nu_nup * d_s1_s2 * G(s1, i) * G(s2, ip) - beta * d_nu_nup * G(s2, m - i - 1) * G(s1, i) + +for s1, s2 in product(spin_names, spin_names): + assert_gfs_are_close(G2_iw_inu_inup_ph_AABB_wick[s1, s2], G2_iw_inu_inup_ph_AABB[s1, s2]) + assert_gfs_are_close(G2_iw_inu_inup_ph_ABBA_wick[s1, s2], G2_iw_inu_inup_ph_ABBA[s1, s2]) + assert_gfs_are_close(G2_iw_inu_inup_pp_AABB_wick[s1, s2], G2_iw_inu_inup_pp_AABB[s1, s2]) + assert_gfs_are_close(G2_iw_inu_inup_pp_ABBA_wick[s1, s2], G2_iw_inu_inup_pp_ABBA[s1, s2])