Developer’s Manual

Programming conventions

Fortran

  1. Indentation is 4 spaces, not tabs. This can be set in your editor.

  2. IMPLICIT NONE is to be used in every Fortran scoping unit.

  3. PRIVATE should be declared within each module, and only those data which need to be available to other modules should be separately declared as PUBLIC.

  4. The USE module, only: syntax of including only what is necessary from other modules should be used to avoid polluting the namespace.

  5. The double precision variables should be declared as REAL(KIND=rKind) for maximum portability.

  6. Generic programming is to be used as much as is reasonable using the strategy outlined in Sec. Generic programming in fortran with templates.

  7. Every loop/if statement in Fortran should be indented.

  8. All Fortran code should be standard conforming and compile using the compilers listed in the manual.

  9. Do not use single-character variable names. Derived types should be capitalized, other variables lower case.

Python

  1. Indentation is 4 spaces, not tabs. This can be set in your editor.

  2. All routines and modules should have a meaningful docstring.

Generic programming in fortran with templates

One of the major drawbacks of programming in Fortran is that it does not have a templating capability, which is the ability of procedures and structures to operate with generic types. Thus, to create an operation scale which scales either real or complex vectors by a real scalar would require the following code stored in a module scale_module.f90:

MODULE scale_module
USE GlobalOps, only: rKind
IMPLICIT NONE

PRIVATE
TYPE Vector
    REAL(KIND=rKind), POINTER :: d(:)
END TYPE Vector

TYPE CVector
    COMPLEX(KIND=rKind), POINTER :: d(:)
END TYPE CVector

 PUBLIC :: scale
 INTERFACE scale
     MODULE PROCEDURE :: scale_Vector, scale_CVector
 END INTERFACE scale

 CONTAINS

 SUBROUTINE scale_Vector(a,v)
     TYPE(Vector) :: v
     REAL(KIND=rKind) :: a

     v%d=a*v%d
 END SUBROUTINE scale_Vector

 SUBROUTINE scale_CVector(a,v)
     TYPE(CVector) :: v
     REAL(KIND=rKind) :: a

     v%d=a*v%d
 END SUBROUTINE scale_CVector
 END MODULE

That is, we essentially have to copy the code twice with only changes to the inputs. This makes the code harder to maintain, and opens up room for type-dependent errors. In MPSFortLib, we have developed a strategy of template programming in Fortran which uses the regex module of Python dealing with regular expressions to build Fortran code from a set of templates. As an example, the module scale_module.f90 above would be replaced by three files. The first would contain a generic template for the vector types, call it scale_types_template.f90:

!!TDic RESET
!!TDic REAL_OR_COMPLEX : real; complex
!!TDic VECTOR_TYPE : vector; cvector

TYPE VECTOR_TYPE
    REAL_OR_COMPLEX(KIND=rKind), POINTER :: d(:)
END TYPE VECTOR_TYPE

Also, a file scale_include.f90 which performs the scaling operation on a VECTOR_TYPE

!!TDic RESET
!!TDic REAL_OR_COMPLEX : real;   complex
!!TDic VECTOR_TYPE :     vector; cvector

 subroutine scale_VECTOR_TYPE(a,v)
     type(VECTOR_TYPE) :: v
     real(KIND=rKind) :: a

     v%d=a*v%d
 end subroutine scale_VECTOR_TYPE
 !PUBLIC_INTERFACE_scale&VECTOR_TYPE

And finally a file scale_template.f90 which instructs how to build scale_module.f90 from these template files:

MODULE scale_module
USE GlobalOps, only: rKind
IMPLICIT NONE

PRIVATE

INCLUDE_TYPES
INCLUDE_PUBLICINTERFACES

CONTAINS

INCLUDE_FILES

END MODULE scale_module

The module scale_module.f90 is built from these three files using our Python code. The Template dictionary (TDic) is defined in front of each type and subroutine after cleaning (!!TDic RESET) it up from previous calls. The word before the colon (:) is the key word to be replaced. The ii-th entry of each line after the color will be used together in the subroutine or function. They are separated by semi-colons. The number of entries after the color specifies the number of subroutines.

The name of a public interface is specified after the subroutine. The name of the interface is before the ampersand &. If the procedure should be public, the statement !PUBLIC_PROCEDURE_scale_VECTOR_TYPE can be added after the subroutine.

This approach allows to use different rules within in template file, which increases the flexibility a bit in comparison to previous version of OSMPS (v1 and v2).

Global keywords always available are

  • !CHECK_ : will disable checks when optimization is choosen.

  • IFMPI_ : enabled when compiling with mpif90. In addition, there are templates for all MPI subroutines.

  • IFIMPS_ and IFNOIMPS_ : will enable (disable) code depending on if iMPS should be compiled or not. Not compiling iMPS drops the dependency on ARPACK.

  • IFMACACC_ and IFNOTMACACC_ : Switch between usage of Mac OS X accelerate package and other operating systems. Right now used to avoid segfaults in BLAS implementation of accelerate package.

For a detailed picture of the regex handling inside Python, please look into …

While all of this machinery is certainly overkill for the scale example above, the nice feature is that the same code scales to an arbitrary number of generic procedures within a module. In the long run, the amount of code and the possibility for errors in duplication is much less.

Building this documentation

This documentation can be build using sphinx, i.e., the standard tool to generate most python documentations. Follow the following steps:

  1. Build the OSMPS library. The documenation will access the fortran modules generated during compilation.

  2. Change to the folder Documenation and run make html

  3. To update the documenation in HTMLDocu, just run a recursive copy cp -r build/html HTMPDocu. This is the folder linked to when building OSMPS.

  4. If you want to update the documentation tracked in the version control, run python3 packdocu.py. Then add and commit to the version control.

Building the local Hilbert space and operators

How MPSPyLib builds Hilbert spaces

To discuss how MPSPyLib builds Hilbert spaces, it is first useful to discuss the notion of lexicographic order. Consider the set \mathcal{V}\left(L,N\right) of vectors with L integer-valued elements with values from 0 to N such that the sum of the elements of the vector is fixed to be N, \sum_{i=1}^Lv_i=N. Take two vectors \mathbf{v} and \mathbf{v}' from this set fulfilling the property that v_i=v_i' if 1\le i\le k-1 and v_k\ne v_k' for some k. We say that \mathbf{v} is superior to \mathbf{v}' if v_k>v_k' and that it is inferior to \mathbf{v}' otherwise. Clearly, \mathbf{v}=\left(N,0,\dots,0\right) is superior to all other vectors and \mathbf{v}=\left(0,\dots,0,N\right) is inferior to all other vectors. Given a vector \mathbf{v}\in\mathcal{V}\left(L,N\right), one finds the next inferior vector in \mathcal{V}\left(L,N\right) by the following algorithm.

  1. If the last element is N then there are no inferior vectors. Exit.

  2. Find the last k such that v_{k+i}=0 for all 1\le i\le L-k.

  3. Set v_k=v_{k}-1, v_{k+1}=N-\sum_{i=1}^{k}v_i, and v_{k+i}=0 for all i>1.

This is the purpose of the function Ops.steplexVec().

Thus, all elements in the set \mathcal{V}\left(L,N\right) may be generated by starting with the completely superior vector \left(N,0,\dots,0\right) and using the above algorithm to generate all other vectors. Other constraints on the minimum or maximum values of the elements of \mathbf{v} are easy to enforce by generating all vectors and checking them individually for the desired properties, and the constraint that the sum be exactly N can be relaxed by using the algorithm for a range of N and adding the spaces in a direct sum. While more efficient algorithms for indexing a Hilbert space may be devised, in MPSPyLib we only build the Hilbert spaces consisting of a single site, and these do not have large dimensionality.

Symmetry-adapted operators

A symmetry-adapted operator is specified as

(10)\hat{O}&=\bigoplus_{\mathbf{q}}\left[\hat{O}_{\mathbf{q}+\Delta \mathbf{q},\mathbf{q}}\right]_{t_{\mathbf{q}+\Delta \mathbf{q}}t_{\mathbf{q}}}\, ,

where \mathbf{q} is the vector of Abelian charges specifying an irrep, t_{\mathbf{q}} is an index running over the degeneracy dimension of irrep \mathbf{q}, and \Delta \mathbf{q} specifies how the operator transforms under the group operation. The definition Eq. (10) assumes that all operators transform covariantly with the group operation. MPSPyLib can not (rather, will not) handle operators which do not transform covariantly. A symmetry-conserving simulation can always be expressed in terms of only covariant operators, and the use of only covariant operators simplifies the numerical implementation and improves its efficiency.

The transformation from operators to their symmetry-adapted counterparts is accomplished with the function named Ops.OperatorList.OperatorstoQOperators(). This function operates with a set of operators expressed as a Python dictionary self.Operators and a matrix or list of generators Generators which are the tags of diagonal matrices whose elements are the charges of the Abelian groups. Since we deal with finite representations, all of the charges can be put into one-to-one correspondence with some subset of the positive integers. This is accomplished by finding a shift \delta for each generator such that the minimum charge is zero and a scaling factor S such that the minimum difference between charges is an integer. That is, we have

G_i&=\left(G_i-\delta_i\right)/S_i

for each generator G_i. The diagonal elements of these new generators define the set of irreps \mathbf{q}, which are sorted to find the associated degeneracy dimensions. All other operators are transformed to this sorted basis. These new operators are stored as a dictionary whose elements are also dictionaries. The keys of the dictionary for a given operator are 'q', a list whose elements are lists of tuples specifying the incoming and outgoing quantum numbers, and 'Blocks', a list of NumPy arrays specifying the action of the operator within the degeneracy space of the specified irrep. In addition, the operator dictionary has as keys 'q_numbers', a list of the irreps of the Hilbert space; 'q_d', a dictionary mapping the irreps to their degeneracy dimensions, 'nqns', the number of Abelian symmetries (i.e. the length of the vector \mathbf{q}); 'q_shifts', the shifts \delta; 'q_scales' the scale factors S; and 'transformation' the matrix transforming the original basis to the symmetry-adapted basis.

Function/Subroutine documentations

The remaining part of the documentation contains the description of all fortran functions and subroutines as far as available. Type definition or interfaces are not considered yet.