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.

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.

Non-template modules

The following modules are not generated by templates:

  • BasicOps
  • GlobalOps
  • ExpokitOps
  • ParallelOps

The other modules are generated before compiling from the set of files located in /templates/.

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 ComplexVector
    COMPLEX(KIND=rKind), POINTER :: d(:)
END TYPE ComplexVector

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

 CONTAINS

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

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

 SUBROUTINE scale_ComplexVector(a,v)
     IMPLICIT NONE
     TYPE(ComplexVector) :: v
     REAL(KIND=rKind) :: a

     v%d=a*v%d
 END SUBROUTINE scale_ComplexVector
 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:

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

scale_VECTOR_TYPE(a,v)
    IMPLICIT NONE
    TYPE(VECTOR_TYPE) :: v
    REAL(KIND=rKind) :: a

    v%d=a*v%d
END SUBROUTINE 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
!PUBLIC_INTERFACE_scale&VECTOR_TYPE

INCLUDE_TYPES
INCLUDE_PUBLICINTERFACES

CONTAINS

INCLUDE_FILES

END MODULE scale_module

The module scale_module.f90 may be built from these three files using the following Python code

INCLUDE_TYPES=''
#Instantiate real types
realdic={
    'VECTOR_TYPE' : 'Vector',
    'REAL_OR_COMPLEX' : 'Real'}

INCLUDE_TYPES+=replace_words_in_File_toString('scale_types_template.f90',realdic)
#Instantiate complex types
cmplxdic={
    'VECTOR_TYPE' : 'ComplexVector',
    'REAL_OR_COMPLEX' : 'Complex'}

INCLUDE_TYPES+=replace_words_in_File_toString('scale_types_template.f90',cmplxdic)
#make a list of the INCLUDE_FILES/Generate them
INCLUDE_FILES=''
INCLUDE_FILES+=replace_words_in_File_toString('scale_include.f90',realdic)
INCLUDE_FILES+=replace_words_in_File_toString('scale_include.f90',cmplxdic)
#Generate the public interfaces
PUBLIC_INTERFACES=ParsePublicInterface([realdic,cmplxdic],'scale_template.f90')
templatedic={
    'INCLUDE_PUBLICINTERFACES' : PUBLIC_INTERFACES,
    'INCLUDE_FILES'     : INCLUDE_FILES,
    'INCLUDE_TYPES'     : INCLUDE_TYPES}
replace_words_in_File_toFile('scale_template.f90','scale_module.f90',templatedic)

First, a dictionary realdic is made which assigns the values Vector and Real to the generic types VECTOR_TYPE and REAL_OR_COMPLEX. This dictionary is then used to replace every instance of VECTOR_TYPE and REAL_OR_COMPLEX in the file scale_types_template.f90 with the corresponding values in realdic and store the result in the string INCLUDE_TYPES. This is accomplished by the routine TemplateGenerator.replace_words_in_File_toString() in MPSPyLib. A similar instantiation is made for the complex vector type using cmplxdic. Next, instantiations of the generic procedures in scale_include.f90 are made using realdic and cmplxdic and stored in the string INCLUDE_FILES. Finally, interfaces to the generic procedures which contain all of the types in realdic and cmplxdic are written to PUBLIC_INTERFACES using the TemplateGenerator.ParsePublicInterface() in MPSPyLib. Finally, the placeholders INCLUDE_PUBLICINTERFACES, INCLUDE_FILES, and INCLUDE_TYPES which appear in scale_template.f90 are replaced with the values INCLUDE\_TYPES, INCLUDE_FILES, and PUBLIC_INTERFACES, respectively, and the resulting code is written to the file scale_module.f90 using the TemplateGenerator.replace_words_in_File_toFile().

The generic routines which are to be combined into a common interface are specified by the syntax

!PUBLIC_INTERFACE_interfaceName&GENERIC_TYPE1&GENERIC_TYPE2&...

Here, interfaceName is the name given to the interface and GENERIC_TYPE1, GENERIC_TYPE2, etc. are the generic types which can be used as arguments to the procedure. In the corresponding *_include.f90 file, the generic functions associated with GENERIC_TYPE1 and GENERIC_TYPE2 should be named interfaceName_GENERIC_TYPE1 and interfaceName_GENERIC_TYPE2, respectively. The code

PUBLIC_INTERFACES=ParsePublicInterface([dic1, dic2,...],'*_template.f90')

will instantiate each generic routine interfaceName according to the instantiations of the associated GENERIC_TYPE in dic1, dic2, etc. A different syntax of the PUBLIC_INTERFACE which is applicable to generic procedures which take multiple generic types as input is

!PUBLIC_INTERFACE_interfaceName&GENERIC_TYPE1,GENERIC_TYPE2&GENERIC_TYPE3,GENERIC_TYPE4&...

Here, GENERIC_TYPE which are separated by a comma indicate multiple generic types which are passed to the same generic procedure. A good example of this functionality is GAXPY in Tensors.f90.

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 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

(1)\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. (1) 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.