# Tag: scientific programming

## Tutorial OOP(V): Documenting Fortran 2003 Classes

In the previous sessions of this tutorial on Object Oriented Programming in Fortran 2003, the basics of OO programming, including the implementation of constructors and destructors as well as operator overloading were covered. The resulting classes have already become quite extended (cf. github source). Although at this point it is still very clear what each part does and why certain choices were made, memory fades. One year from now, when you revisit your work, this will no longer be the case. Alternately, when sharing code, you don’t want to have to dig through every line of code to figure out how to use it. These are just some of the reasons why code documentation is important. This is a universal habit of programming which should be adopted irrespective of the programming-language and-paradigm, or size of the code base (yes, even small functions should be documented).

In Fortran, comments can be included in a very simple fashion: everything following the “!” symbol (when not used in a string) is considered a comment, and thus ignored by the compiler. This allows for quick and easy documentation of your code, and can be sufficient for single functions. However, when dealing with larger projects retaining a global overview and keeping track of interdependencies becomes harder. This is where automatic documentation generation software comes into play.  These tools parse specifically formatted comments to construct API documentation and user-guides. Over the years, several useful tools have been developed for the Fortran language directly, or as a plugin/extension to a more general tool:

• ROBODoc : A tool capable of generating documentation (many different formats) for any programming/script language which has comments.  The latest update dates from 2015.
• Doctran : This tool is specifically aimed at free-format (≥ .f90 ) fortran, and notes explicitly the aim to deal with object oriented f2003. It only generates html documentation, and is currently proprietary with license costs of 30£ per plugin. Latest update 2016.
• SphinxFortran : This extension to SphinxFortran generates automatic documentation for f90 source (no OO fortran) and generates an html manual. This package is written in python and requires you to construct your config file in python as well.
• f90doc / f90tohtml : Two tools written in Perl, which transform f90 code into html webpages.
• FotranDOC : This tool (written in Fortran itself) aims to generate documentation for f95 code, preferably in a single file, in latex. It has a simple GUI interface, and the source of the tool itself is an example of how the fortran code should be documented. How nice is that?
• FORD : Ford is a documentation tool written in python, aimed at modern fortran (i.e. ≥ f90).
• Doxygen :  A multi-platform automatic documentation tool developed for C++, but extended to many other languages including fortran. It is very flexible, and easy to use and can produce documentation in html, pdf, man-pages, rtf,… out of the box.

As you can see, there is a lot to choose from, all with their own quirks and features. One unfortunate aspect is the fact that most of these tools use different formatting conventions, so switching from one to the another is not an exercise to perform lightly. In this tutorial, the doxygen tool is used, as it provides a wide range of options, is multi-platform,  supports multiple languages and multiple output formats.

As you might already expect, Object Oriented Fortran (f2003) is a bit more complicated to document than  procedural Fortran, but with some ingenuity doxygen can be made to provide nice documentation even in this case.

## 1. Configuring Doxygen

Before you can start you will need to install doxygen:

1. Go the the doxygen-download page and find the distribution which is right for you (Windows-users: there are binary installers, no hassle with compilations 🙂 ).
2. Follow the installation instructions, also install GraphViz, this will allow you to create nicer graphics using the dot-tool.
3. Also get a pdf version of the manual (doxygen has a huge number of options)

With a nicely installed doxygen, you can make use of the GUI to setup a configuration suited to your specific needs and generate the documentation for your code automatically. For Object Oriented Fortran there are some specific settings you should consider:

1. #### Wizard tab

• Project Topic : Fill out the different fields. In a multi-file project, with source stored in a folder structure, don’t forget to select the tick-box “Scan recursively” .
• Mode Topic : Select “Optimize for Fortran output”.
• Output Topic : Select one or more output formats you wish to generate: html, Latex (pdf), map-pages, RTF, and XML
• Diagrams Topic: Select which types of diagrams you want to generate.
2. #### Expert tab

(Provides access each single configuration option to set in doxygen, so I will only highlight a few. Look through them to get a better idea of the capabilities of doxygen.)

• Project Topic :
• EXTENSION_MAPPING: You will have to tell doxygen which fortran extensions you are using by adding them, and identifying it as free format fortran: e.g. f03=FortranFree (If you are also including text-files to provide additional documentation, it is best to add them here as well as free format fortran).
• Build Topic:
• CASE_SENSE_NAMES: Even though Fortran itself is not case sensitive, it may be nice to keep the type of casing you use in your code in your documentation. Note, however, that even though the output may have upper-case names, the documentation itself will require lower-case names in references.
• Messages Topic:
• WARN_NO_PARAMDOC: Throw a warning if documentation is missing for a function variable. This is useful to make sure you have a complete documentation.
• Source Browser Topic:
• SOURCE_BROWSER: Complete source files are included in the documentation.
• INLINE_SOURCES: Place the source body with each function directly in the documentation.
• HTML Topic:
• FORMULA_FONTSIZE: The fontsize used for generated formulas. If 10 pts is too small to get a nice effect of formulas embedded in text.
• Dot Topic:
• HAVE_DOT & DOT_PATH: If you installed GraphViz
• DOT_GRAPH_MAX_NODES: Maximum number of nodes to draw in a relation graph. In case of larger projects, 50 may be too small.
• CALL_GRAPH & CALLER_GRAPH: Types of relation graphs to include.
3. #### Run tab

• Press “Run doxygen” and watch how your documentation is being generated. For larger projects this may take some time. Fortunately, graphics are not generated anew if they are present from a previous run, speeding things up. (NOTE: If you want to generate new graphics (and equations with larger font size), make sure to delete the old versions first.) Any warnings and errors are also shown in the main window.
• Once doxygen was run successfully, pressing the button “Show HTML output” will open a browser and take you to the HTML version of the documentation.

Once you have a working configuration for doxygen, you can save this for later use. Doxygen allows you to load an old configuration file and run immediately. The configuration file for the Timer-class project is included in the docs folder, together with the pdf-latex version of the generated documentation.  Doxygen generates all latex files required for generating the pdf. To generate the actual pdf, a make.bat file needs to be run (i.e. double-click the file, and watch it run) in a Windows environment.

## 2. Documenting Fortran (procedural)

Let us start with some basics for documenting Fortran code in a way suitable for doxygen. Since doxygen has a very extensive set of options and features, not all of them can be covered. However, the manual of more than 300 pages provides all the information you may need.

With doxygen, you are able to document more or less any part of your code: entire files, modules, functions or variables. In each case, a similar approach can be taken. Let’s consider the documentation of the TimeClass module:

 Documentation of the TimeClass module.
1. !++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
2. !> \brief The <b>TimeClass module</b> contains the
5. !!
6. !! @author  Dr. Dr. Danny E. P. Vanpoucke
7. !! @version 2.0-3  (upgrades deprecated timing module)
8. !! @date    19-03-2020
10. !!
11. !! @warning Internally, Julian Day Numbers are used to compare dates. As a
12. !! result, *negative* dates are not accepted. If such dates are created
13. !! (*e.g.*, due to a subtraction), then the date is set to zero.
14. !!
15. !! This module makes use of:
16. !! - nothing; this module is fully independent
17. !<-----------------------------------------------------------------------
18. module TimeClass
19.     implicit none
20.     private

The documentation is placed in a standard single or multi-line fortran comment.  In case of multi-line documentation, I have the personal habit turning it into a kind of banner starting with a “!+++++++++” line and closing with a “!<——————-” line.  Such choices are your own, and are not necessary for doxygen documentation. For doxygen, a multi-line documentation block starts with “!>” and ends with “!<“ . The documentation lines in between can be indicated with “!!”. This is specifically for fortran documentation in doxygen. C/C++ and other languages will have slightly different conventions, related to their comment section conventions.

In the block above, you immediately see certain words are preceded by an “@”-symbol or a “\”, this indicates these are special keywords. Both the “@” and “\” can be used interchangeably for most keywords, the preference is again personal taste.  Furthermore, doxygen supports both html and markdown notation for formatting, providing a lot of flexibility. The multi-line documentation is placed before the object being documented (here an entire module).

Some keywords:

• \brief : Here you can place a short description of the object. This description is shown in parts of the documentation that  provide an overview. Note that this is also the first part of the full documentation of the object itself. After a blank line, the \details(this keyword does not need to provided explicitly) section starts, providing further details on the object. This information is only visible in the documentation of the object itself.
• “::”  : Referring to an element of an object can be done by linking the element and the object via two colons:  object::element . Here it is important to remember that your module is an object, so linking to an element of a module from outside that module requires you to refer to it in this way.
• @author : Provide information on the author.
• @version : Provide version information.
• @date : Provide information on the date.
• @warning : Provides a highlighted section with warning information for the user of your code (e.g., function kills the program when something goes wrong).
• @todo[not shown] If you still have some things to do with regard to this object you can use this keyword. More interestingly, doxygen will also create a page where all to-do’s of the entire project are gathered, and link back to the specific code fragments.

 Documentation of a function
1. !++++++++++++++++++++++++++++++++++++++++++++++
3. !! via the "-" operator. This is the function
5. !!
6. !! \b usage:
7. !! \code{.f03}
8. !! Total = this - that
9. !! \endcode
11. !!
12. !! \note The result should remain a positive number.
13. !!
15. !!                 the "-" operator.
17. !!                 the "-" operator.
19. !!               the difference.
20. !<---------------------------------------------
21.     pure function subtract(this,that) Result(Total)
22.         class(TTime), intent(in) :: this, that
23.         Type(TTime) :: total

When documenting functions and subroutines there are some addition must-have keywords.

• @param[in] , @param[out] ,or@param[in,out] : Provide a description for each of the function parameters, including their  intent: “in”, “out”, or “in,out” (note the comma!).
• \return : Provides information on the return value of the function.
• \b, \i : The next word is bold or italic
• \n : Start a newline, without starting a new paragraph.
• \note : Add a special note in your documentation. This section will be high lighted in a fashion similar to @warning.
• \code{.f03}…\endcode :  This environment allows you to have syntax highlighted code in your documentation. The language can be indicated via the “extension” typical for said language. In this case: fortran-2003.
• \f$… \f$, or \f[ … \f] : Sometimes equations are just that much easier to convey your message. Doxygen also supports latex formatting for equations. These tags can be used to enter a latex $…$ or  math environments. The equations are transformed into small png images upon documentation generation, to be included in the html of your documentation. There are two important aspects to consider when using this option:
1. Font size of the equation: Check if this is sufficient and don’t be afraid to change the font size to improve readability.
2. Compilation is not halted upon an error: If the latex compiler encounters an error in your formula it just tries to continue. In case of failure, the end result may be missing or wrong. Debugging latex equations in doxygen documentation can be quite challenging as a result. So if you are using large complex equations, it may be advised to run them in a pure latex environment, and only past them in the documentation once you are satisfied with the result.

## 3. Documenting Fortran Classes

With the knowledge of the previous section, it is relatively easy to document most fortran code. Also the type of object orientation available in fortran 95, in which a fortran module is refurbished as a class. True fortran classes in contrast tend to give a few unexpected issues to deal with. Lets have a look at the documentation of the TTime class of the TimeClass module:

 Documentation of a fortran class definition
1. !+++++++++++++++++++++++++++++++++++++++
2. !> @class ttime
3. !! \brief The TTime class contains all time functionality
4. !! with regard to a single time stamp.
5. !<-------------------------------------
6.     type, public :: TTime
7.       private
8.         integer :: year    !< @private The year
9.         integer :: month   !< @private The month (as integer).
10.         ...
11.     contains
12.       private
13.         procedure, pass(this),public :: SetTime       !<          @copydoc timeclass::settime
14.         procedure, pass(this)        :: CalculateJDN  !< @private @copydoc timeclass::calculatejdn
15.         procedure, pass(this)        :: SetJDN        !< @private @copydoc timeclass::setjdn
16.         ...
17.         procedure, pass(this)        :: copy          !< @private @copydoc timeclass::copy
18.         ...
19.         generic, public :: assignment(=) => copy      !<          @copydoc timeclass::copy
20.         !> @{ @protected
21.         final :: destructor !< @copydoc timeclass::destructor
22.         !> @}
23.     end type TTime
24.
25.     ! This is the only way a constructor can be created,
26.     ! as no "initial" exists, emulates the C++ constructor behavior
27.     interface TTime
28.         module procedure constructor
29.     end interface TTime

To make sure doxygen generates a class-like documentation for our fortran class, it needs to be told it is a class. This can be done by documenting the class itself and using the keyword @class nameclass, with nameclass the name doxygen will use for this class (so you can choose something different from the actual class name). Unfortunately, doxygen will call this a “module” in the documentation (just poor luck in nomenclature). On the module page for the ttime class a listing is provided of all elements given in the class definition. The documentation added to each member (e.g.,:

 Source code
1. integer :: year !< @private The year

is shown as “\brief” documentation. By default all members of our function are considered as public. Adding the @private, @public, or @protected keyword instructs doxygen explicitly to consider these members as private, public or protected. (I used protected in the ttime code not as it should be used in fortran, but as a means of indicating the special status of the final subroutine (i.e. protected in a C++ way).)

However, there seems to be something strange going on. When following the links in the documentation, we do not end up with the documentation provided for the functions/subroutines in the body of our timeclass module. Doxygen seems to consider these two distinct things. The easiest way to link the correct information is by using the keyword @copydoc functionreference . The documentation is (according to doxygen) still for two distinctly different objects, however, this time they have the exact same documentation (unless you add more text on the member documentation line). In this context, it interesting to know there is also @copybrief and @copydetails which can be used to only copy the brief/details section.

In this example, the constructor interface is not documented, as this created confusion in the final  documentation since doxygen created a second ttime module/object linked to this interface. However, not documenting this specific instance of the constructor does not create such a large issue, as the module(the fortran module) function itself is documented already.

## Conclusion

Documenting fortran classes can be done quite nicely with doxygen. It provides various modes of output: from a fully working website with in-site search engine to a hyperlinked pdf or RTF document. The flexibility and large number of options may be a bit daunting at first, but you can start simple, and work your way up.

As Fortran is supported as an extension, you will need to play around with the various options to find which combination gives the effect you intended. This is an aspect present in all automated code documentation generation tools, since object oriented Fortran is not that widely used. Nonetheless, doxygen provides a very powerful tool worth your time and effort.

#### PS:

In the previous tutorial, we created a constructor and destructor for our TTimer class.  Next, we extend our class with overloaded operators. Depending on the type of object your class represents, you may want to define an addition/subtraction/multiplication/… operator. In addition, the assignment operator deserves some extra attention as you may want to have a clear control over this operation  (e.g.deep copy vs shallow copy). The full source of this tutorial and the previous, can be downloaded from my github-page.

Let us start with the latter: the assignment operator. As with all other operators, it is possible to overload the assignment operator in modern fortran.

When dealing with objects and classes—or extended data-structures in general—, their properties often are (implicit) pointers to the actual data-structure. This brings an interesting source of possible bugs due to shallow copies being made while deep copies are expected (although the problem may be less pronounced in Fortran than it is in Python).

In a fortran object, the assignment of a pointer component (i.e., an explicit pointer variable, or a component which is an object itself) happens via a shallow copy (or pointer assignment). In contrast, for an allocatable component, the assignment operation performs by default a deep copy (i.e., space is allocated, and values are copied). Shallow copies are very useful with regard to quickly creating new handles to the same data-structure. However, if you want to make a true copy, which you can modify without changing the original, then a deep copy is what you want. By implementing assignment overloading for your own classes, you have more control over the actual copying process, and you can make sure you are creating deep copies if those are preferred.

The implementation of overloading for the assignment operator is not too complicated. It requires two lines in your class definition:

type, public :: TTimer
private
...
contains
private
procedure, pass(this) :: Copy                   !< Make a copy of a timer object
generic, public       :: assignment(=) => Copy  !< This is how copy is used.
...
end type TTimer

First, you need to define a class method which performs a copy-operation—which in a fit or original though we decided to call “copy” ;-).  As you can see this function is private, so it will not be accessible to the user of your class via a call like :

call MyTimer%Copy()

Secondly, you link this class method via the “=>” to the assignment-operator.  It is a generic interface, which means the assignment operator could be linked to different functions, of which the relevant one will be determined and used during run-time. This generic is also public  (otherwise you would not be able to use it).

The implementation of the class method follows the standard rules of any class method and could look like

pure subroutine Copy(this,from)
class(TTimer), intent(inout) :: this
class(TTimer), intent(in) :: from

this%firstProperty = from%firstProperty
...
!make explicit copies of all properties and components
...

end subroutine Copy

The “this” object which we passed to our class method is the object on the left side of the assignment operator, while the “from” object is the one on the right side. Note that both objects are defined as “class” and not as “type”. Within the body of this method you are in charge of copying the data from the “from”-object to the “this”-object, giving you control over deep/shallow copying.

In practice the overloaded operator is used as:

type(TTimer):: TimerThis, TimerFrom

TimerFrom = TTimer() ! initialization of the timers
TimerThis = TTimer() ! (cf., previous tutorial on constructors and destructors)
...
! do stuff with TimerFrom
...
TimerThis = TimerFrom ! although you type "=", the overloading causes this to be implemented as-if you wrote
! call TimerThis%copy(TimerFrom)

Just as you can overload the assignment operator above, you can also overload all other fortran operators. However, be careful to keep things intuitive.  For example, an addition operation on our TTimer class is strange. What would it mean to add one timer to another? How would you subtract one chronometer from another? In contrast, inside our TTimer class we have a list of TTime objects which can be used to represent a date and time, as-well-as a time interval.[1]  For the remainder of this tutorial, we will assume the TTime class only represents time-intervals. For such a class, it makes sense to be able to add and subtract time intervals.

type, public :: TTime
private
...
! the properties of the TTime class
...
contains
private
...
! the methods of the TTime class
...
procedure, pass(this)        :: copy          ! Copy content from other TTime instance,
! private, accessed via the assignment statement
procedure, pass(this)        :: subtract      ! subtract two TTime instances.
generic, public :: assignment(=) => copy      ! This is how copy is used.
generic, public :: operator(+)   => add       ! This is how add is used.
generic, public :: operator(-)   => subtract  ! This is how subtract is used.
final :: destructor
end type TTime

interface TTime
module procedure constructor
end interface TTime


pure function add(this,that) Result(Total)
class(TTime), intent(in) :: this, that
Type(TTime) :: total

total = TTime()
...
! implementation of the addition of the properties of
! this to the properties of that, and storing them in
! Total
! e.g.: Total%seconds = this%seconds + that%seconds
...


The returned object need to be defined as a type, and the further implementation of the function follows the standard fortran rules. It is important to note that for a function-header like this one, the object to the left of the operator will be the one calling the overloaded operator function, so:

Total = this + that

and not

Total = that + this

This may not seem this important, as we are adding two objects of the same class, but that is not necessarily always the case. Imagine that you want to overload the multiplication operator, such that you could multiply your time-interval with any possible real value. On paper

Δt * 3.5 = 3.5 * Δt

but for the compiler in the left product “this” would be a TTime object and “that” would be a real, while in the right product “this” is the real, and “that” is the TTime object. To deal with such a situation, you need to implement two class methods, which in practice only differ in their header:

pure function MultLeft(this,that) Result(Total)
class(TTime), intent(in) :: this
real, intent(in) :: that
Type(TTime) :: total

and

pure function MultRight(that, this) Result(Total)
class(TTime), intent(in) :: this
real, intent(in) :: that
Type(TTime) :: total

In the class definition both functions are linked to the operator as

procedure, pass(this) ::  MultLeft
procedure, pass(this) ::  MultRight
generic, public :: operator(*) => MultLeft, MultRight

With this in mind, we could also expand our implementation of the “+” and “” operator, by adding functionality that allows for the addition and subtraction of reals representing time-intervals. Also here, the left and right versions would need to be implemented.

As you can see, modern object oriented fortran provides you all the tools you need to create powerful classes capable of operator overloading using simple and straightforward implementations.

In our next Tutorial, we’ll look into data-hiding and private/public options in fortran classes.

[1] You could argue that this is not an ideal choice and that it would be better to keep these two concepts ( absolute and relative time) separate through the use of different classes.

## Tutorial OOP(III): Constructors and Destructors

In this tutorial on Object Oriented Programming in Fortran 2003, we are going to discuss how to create constructors and destructors for a Fortran class. During this tutorial, I assume that you know how to create a new project and what a class looks like in Fortran 2003.  This tutorial is build around a TimerClass, which I wrote as an upgrade for my initial timing module in HIVE-tools. The full source of this TimerClass can be found and downloaded from github.

Where the former two tutorials were aimed at translating a scientific model into classes within the confines of the Fortran programming language, this tutorial is aimed at consolidating a class using good practices: The creation a constructor and destructor. As the destructor in Fortran classes is most straight forward of the two, we’ll start with it.

## 1. The destructor.

A destructor is a method (i.e., a class subroutine) which is automatically invoked when the object is destroyed (e.g., by going out of scope).  In case of a Fortran class, this task is performed by the class-method(s) indicated as  final procedure. Hence such methods are also sometimes referred to as finalizers. Although in some languages destructors and finalizers are two distinctly different features (finalizers are then often linked to garbage collecting), within the Fortran context I consider them the same.

Within the definition of our TTimerClass the destructor is implemented as:

 Destructor of the TTimerClass
1. module TimerClass
2. implicit none
3.
4.     type, public :: TTimer
5.       private
6.       ! here come the properties
7.     contains
8.       private
9.       ! here come the methods
10.       final :: destructor
11.     end type TTimer
12.
13. contains
14.
15.     subroutine destructor(this)
16.     Type(TTimer) :: this
17.     ! Do whatever needs doing in the destructor
18.     end subroutine destructor
19.
20. end module TimerClass

In contrast to a normal class-method, the destructor is called using the final keyword, instead of the usual procedure keyword. This method is private, as it is not intended to be used by the user anyway, only by the compiler upon cleanup of the instance of the class (i.e., the object). Furthermore, although defined as part of the class, a final subroutine is not type-bound, and can thus not be accessed through the type.

The destructor subroutine itself is a normal Fortran subroutine. There is, however, one small difference with a usual class-method, the parameter referring to the object (c.q. “this“) is indicated as a TYPE and not as a CLASS. This is because the destructor is only applicable to properties belonging to this “class” (Note that final subroutines are not inherited by the child-class). For a child-class (also called a derived class), the destructor of the child-class should deal with all the additional properties of the child-class, while the destructor of the parent-class is called to deal with its respective properties. In practice, the destructor of the child-class is called first, after which the destructor of the parent class is called (and recursively further along the class its family tree.)

So what do you put in such a destructor? Anything that needs to be done to allow the object to be gracefully terminated. Most obviously: deallocation of allocatable arrays, pointer components, closing file handles,…

## 2. The constructor.

Where other programming  languages may provide an initialization section or access to a key-worded constructor. Although Fortran allows for variables to be initialized upon definition, there is no constructor keyword available to be used in its classes. Of course, this does not prevent you from adding an “init()” subroutine which the user should call once the new object is allocated. You could even use a private Boolean property (initialized old style)  to keep track of the fact that an object was initialized when entering any of its methods, and if not, call the init() function there and then. There are many ways to deal with the initialization of a new object.  Furthermore, different approaches also put the burden of doing things right either with the programmer developing the class, or the user, applying the class and creating objects.

Here, I want to present an approach which allows you to present a clear set-up of your class and which resembles the instance creation approach also seen in other languages (and which implicitly shows the “pointer”-nature of objects ):

NewObject = TClass()

In case of our TTimer class this will look like:

Type(TTimer) :: MyTimer

MyTimer = TTimer()

This means we need to have a function with the exact same name as our class (cf., above), which is achieved through the use of an interface to a module procedure.  Just giving this name to the constructor function itself will cause your compiler to complain (“Name ttimer at (1)  is already defined as a generic interface“).  By using a different name for the  function, and wrapping it in an interface, this issue is avoided.

 Class Constructor
1. module Timerclass
2.     implicit none
3.
4.
5.     type, public :: TTimer
6.         private
7.     ...
8.     contains
9.         private
10.         ...
11.     end type TTimer
12.
13.     interface TTimer
14.         module procedure Constructor
15.     end interface TTimer
16.
17. contains
18. function Constructor() Result(Timer)
19.     type(TTimer) :: Timer
20.
21.     !initialize variables directly
22.     Timer%x=...
23.     ! or through method calls
24.     call Timer%setTime(now)
25.     ...
26.
27. end function Constructor
28.
29. end module TimerClass

Note that the constructor function is not part of the class definition, and as such the object is not passed to the constructor function. In addition, the Timer object being created is defined as a Type(TTimer) not Class(TTimer), also because this function is not part of the class definition.

That is all there is to it. Simple and elegant.

In our next Tutorial, we’ll have a look at operator and assignment overloading. Combined with a constructor and destructor as presented here, you are able to create powerful and intuitive classes (even in Fortran).

## New year’s resolution

A new year, a new beginning.

For most people this is a time of making promises, starting new habits or stopping old ones. In general, I forgo making such promises, as I know they turn out idle in a mere few weeks without external stimulus or any real driving force.

In spite of this, I do have a new years resolution for this year: I am going to study machine learning and use it for any suitable application I can get my hands on (which will mainly be materials science, but one never knows).  I already have a few projects in mind, which should help me stay focused and on track. With some luck, you will be reading about them here on this blog. With some more luck, they may even end up being part of an actual scientific publication.

But first things first, learn the basics (beyond hear-say messages of how excellent and world improving AI is/will be). What are the different types of machine learning available, is it all black box or do you actually have some control over things. Is it a kind of magic? What’s up with all these frameworks (isn’t there anyone left who can program?), and why the devil seem they all to be written in a script langue (python) instead of a proper programming language? A lot of questions I hope to see answered. A lot of things to learn. Lets start by building some foundations…the old fashioned way: By studying using a book, with real paper pages!

## Defensive Programming and Debugging

Last few months, I finally was able to remove something which had been lingering on my to-do list for a very long time: studying debugging in Fortran. Although I have been programming in Fortran for over a decade, and getting quite good at it, especially in the more exotic aspects such as OO-programming, I never got around to learning how to use decent debugging tools. The fact that I am using Fortran was the main contributing factor. Unlike other languages, everything you want to do in Fortran beyond number-crunching in procedural code has very little documentation (e.g., easy dll’s for objects), is not natively supported (e.g., find a good IDE for fortran, which also supports modern aspects like OO, there are only very few who attempt), or you are just the first to try it (e.g., fortran programs for android :o, definitely on my to-do list). In a long bygone past I did some debugging in Delphi (for my STM program) as the debugger was nicely integrated in the IDE. However, for Fortran I started programming without an IDE and as such did my initial debugging with well placed write statements. And I am a bit ashamed to say, I’m still doing it this way, because it can be rather efficient for a large code spread over dozens of files with hundreds of procedures.

However, I am trying to repent for my sins. A central point in this penance was enlisting for the online MOOC  “Defensive Programming and Debugging“. Five weeks of intense study followed, in which I was forced to use command line gdb and valgrind. During these five weeks I also sharpened my skills at identifying possible sources of bugs (and found some unintentional bugs in the course…but that is just me). Five weeks of hard study, and taking tests, I successfully finished the course, earning my certificate as defensive programmer and debugger. (In contrast to my sometimes offensive programming and debugging skills before 😉 .)

## A Spectre and Meltdown victim: VASP

Over the last weekend, two serious cyber security issues were hot news: Meltdown and Spectre [more links, and links](not to be mistaken for a title of a bond-movie). As a result, also academic HPC centers went into overdrive installing patches as fast as possible. The news of the two security issues went hand-in-hand with quite a few belittling comments toward the chip-designers ignoring the fact that no-one (including those complaining now) discovered the problem for over decade. Of course there was also the usual scare-mongering (cyber-criminals will hack our devices by next Monday, because hacks using these bugs are now immediately becoming their default tools etc.) typical since the beginning of the  21st century…but now it is time to return back to reality.

One of the big users on scientific HPC installations is the VASP program(an example), aimed at the quantum mechanical simulation of materials, and a program central to my own work. Due to an serendipitous coincidence of a annoyingly hard to converge job I had the opportunity to see the impact of the Meltdown and Spectre patches on the performance of VASP: 16% performance loss (within the range of the expected 10-50% performance loss for high performance applications [1][2][3]).

The case:

• large HSE06 calculation of a 71 atom defective ZnO supercell.
• 14 irreducible k-points (no reduction of the Hartree-Fock k-points)
• 14 nodes of 24 cores, with KPAR=14, and NPAR=1 (I know NPAR=24 is the recommended option)

The calculation took several runs of about 10 electronic steps (of each about 5-6 h wall-time, about 2.54 years of CPU-time per run) . The relative average time is shown below (error-bars are the standard deviation of the times within a single run). As the final step takes about 50% longer it is treated separately. As you can see, the variation in time between different electronic steps is rather small (even running on a different cluster only changes the time by a few %). The impact of the Meltdown/Spectre patch gives a significant impact.

Impact of Meltdown/Spectre patch on VASP performance for a 336 core MPI job.

The HPC team is currently looking into possible workarounds that could (partially) alleviate the problem. VASP itself is rather little I/O intensive, and a first check by the HPC team points toward MPI (the parallelisation framework required for multi-node jobs) being ‘a’ if not ‘the’ culprit. This means that also an impact on other multi-node programs is to be expected. On the bright side, finding a workaround for MPI would be beneficial for all of them as well.

So far, tests I performed with the HPC team not shown any improvements (recompiling VASP didn’t help, nor an MPI related fix). Let’s keep our fingers crossed, and hope the future brings insight and a solution.

## Exa-scale computing future in Europe?

As a computational materials scientist with a main research interest in the ab initio simulation of materials, computational resources are the life-blood of my research. Over the last decade, I have seen my resource usage grow from less than 100.000 CPU hours per year to several million CPU-hours per year. To satisfy this need for computational resources I have to make use of HPC facilities, like the TIER-2 machines available at the Flemish universities and the Flemish TIER-1 supercomputer, currently hosted at KU Leuven. At the international level, computational scientists have access to so called TIER-0 machines, something I no doubt will make use of in the future. Before I continue, let me first explain a little what this TIER-X business actually means.

The TIER-X notation is used to give an indication of the size of the computer/supercomputer indicated. There are 4 sizes:

•  TIER-3: This is your personal computer(laptop/desktop) or even a small local cluster of a research group. It can contain from one (desktop) up to a few hundred CPU’s (local cluster). Within materials research, this is sufficient for quite a few tasks: post-processing of data, simple force-field based calculations, or even small quantum chemical or solid state calculations. A large fraction of the work during my first Ph.D. was performed on the local cluster of the CMS.
• TIER-2: This is a supercomputer hosted by an institute or university. It generally contains over 1000 CPUs and has a peak performance of >10 TFLOPS (1012 Floating Point Operations Per Second, compare this to 1-50×10FLOPS or 1-25 GFLOPS of an average personal computer). The TIER-2 facilities of the VUB and UAntwerp both have a peak performance of about 75TFLOPS , while the machines at Ghent University and the KU Leuven/Uhasselt facilities both have a peak performance of about 230 TFLOPS. Using these machines I was able to perform the calculations necessary for my study of dopant elements in cerates (and obtain my second Ph.D.).
• TIER-1: Moving up one more step, there are the national/regional supercomputers. These generally contain over 10.000 CPUs and have a peak performance of over 100 TFLOPS. In Flanders the Flemish Supercomputer Center (VSC) manages the TIER-1 machine (which is being funded by the 5 Flemish universities). The first TIER-1 machine was hosted at Ghent University, while the second and current one is hosted at KU Leuven, an has a peak performance of 623 TFLOPS (more than all TIER-2 machines combined), and cost about 5.5 Million € (one of the reasons it is a regional machine). Over the last 5 years, I was granted over 10 Million hours of CPU time, sufficient for my study of Metal-Organic Frameworks and defects in diamond.
• TIER-0: This are international level supercomputers. These machines contain over 100.000 CPUs, and have a peak performance in excess of 1 PFLOP (1 PetaFLOP = 1000 TFLOPS). In Europe the TIER-0 facilities are available to researchers via the PRACE network (access to 7 TIER-0 machines, accumulated 43.49 PFLOPS).

This is roughly the status of what is available today for Flemish scientists at various levels. With the constantly growing demand for more processing power, the European union, in name of EuroHPC, has decided in march of this year, that Europe will host two Exa-scale computers. These machines will have a peak performance of at least 1 EFLOPS, or 1000 PFLOPS. These machines are expected to be build by 2024-2025. In June, Belgium signed up to EuroHPC as the eighth country participating, in addition, to the initial 7 countries (Germany, France, Spain, Portugal, Italy, Luxemburg and The Netherlands).

This is very good news for all involved in computational research in Flanders. There is the plan to build these machines, there is a deadline, …there just isn’t an idea of what these machines should look like (except: they will be big, massively power consuming and have a target peak performance). To get an idea what users expect of such a machine, Tier-1 and HPC users have been asked to put forward requests/suggestions of what they want.

From my user personal experience, and extrapolating from my own usage I see myself easily using 20 million hours of CPU time each year by the time these Exa-scale machines are build. Leading a computational group would multiply this value. And then we are talking about simple production purpose calculations for “standard” problems.

The claim that an Exa-scale scale machine runs 1000x faster than a peta-scale machine, is not entirely justified, at least not for the software I am generally encountering. As software seldom scales linearly, the speed-gain from Exa-scale machinery mainly comes from the ability to perform many more calculations in parallel. (There are some exceptions which will gain within the single job area, but this type of jobs is limited.) Within my own field, quantum mechanical calculation of the electronic structure of periodic atomic systems, the all required resources tend to grow with growth of the problem size. As such, a larger system (=more atoms) requires more CPU-time, but also more memory. This means that compute nodes with many cores are welcome and desired, but these cores need the associated memory. Doubling the cores would require the memory on a node to be doubled as well. Communication between the nodes should be fast as well, as this will be the main limiting factor on the scaling performance. If all this is implemented well, then the time to solution of a project (not a single calculation) will improve significantly with access to Exa-scale resources. The factor will not be 100x from a Pflops system, but could be much better than 10x. This factor 10 also takes into account that projects will have access to much more demanding calculations as a default (Hybrid functional structure optimization instead of simple density functional theory structure optimization, which is ~1000x cheaper for plane wave methods but is less accurate).

At this scale, parallelism is very important, and implementing this into a program is far from a trivial task. As most physicists/mathematicians/chemists/engineers may have the skills for writing scientifically sound software, we are not computer-scientists and our available time and skills are limited in this regard. For this reason, it will become more important for the HPC-facility to provide parallelization of software as a service. I.e. have a group of highly skilled computer scientists available to assist or even perform this task.

Next to having the best implementation of software available, it should also be possible to get access to these machines. This should not be limited to a happy few through a peer review process which just wastes human research potential. Instead access to these should be a mix of guaranteed access and peer review.

• Guaranteed access: For standard production projects (5-25 million CPU hours/year) university researchers should have a guaranteed access model. This would allow them to perform state of the are research without too much overhead. To prevent access to people without the proven necessary need/skills it could be implemented that a user-database is created and appended upon each application. Upon first application, a local HPC-team (country/region/university Tier-1 infrastructure) would have to provide a recommendation with regard to the user, including a statement of the applicant’s resource usage at that facility. Getting resources in a guaranteed access project would also require a limited project proposal (max 2 pages, including user credentials, requested resources, and small description of the project)
• Peer review access: This would be for special projects, in which the researcher requires a huge chunk of resources to perform highly specialized calculations or large High-throughput exercises (order of 250-1000 million CPU hours, e.g. Nature Communications 8, 15959 (2017)). In this case a full project with serious peer review (including rebuttal stage, or the possibility to resubmit after considering the indicated problems). The goal of this peer review system should not be to limit the number of accepted projects, but to make sure the accepted projects run successfully.
• Pay per use: This should be the option for industrial/commercial users.

What could an HPC user as myself do to contribute to the success of EuroHPC? This is rather simple, run the machine as a pilot user (I have experience on most of the TIER-2 clusters of Ghent University and both Flemish Tier-1 machines. I successfully crashed the programs I am using by pushing them beyond their limits during pilot testing, and ran into rather unfortunate issues. 🙂 That is the job of a pilot user, use the machine/software in unexpected ways, such that this can be resolved/fixed by the time the bulk of the users get access.) and perform peer review of the lager specialized projects.

Now the only thing left to do is wait. Wait for the Exa-scale supercomputers to be build…7 years to go…about 92 node-days on Breniac…a starting grant…one long weekend of calculations.

##### Appendix

For simplicity I use the term CPU to indicate a single compute core, even though technically, nowadays a single CPU will contain multiple cores (desktop/laptop: 2-8 cores, HPC-compute node: 2-20 cores / CPU (or more) ). This to make comparison a bit more easy.

Furthermore, modern computer systems start more and more to rely on GPU performance as well, which is also a possible road toward Exa-scale computing.

Orders of magnitude:

• G = Giga = 109
• T = Tera = 1012
• P = Peta = 1015
• E = Exa = 1018

## Resource management on HPC infrastructures.

Computational as a third pillar of science (next to experimental and theoretical) is steadily developing in many fields of science. Even some fields you would less expect it, such as sociology or psychology. In other fields such as physics, chemistry or biology it is much more widespread, with people pushing the boundaries of what is possible. Larger facilities provide access to larger problems to tackle. If a computational physicist is asked if larger infrastructures would not become too big, he’ll just shrug and reply: “Don’t worry, we will easily fill it up, even a machine 1000x larger than that.” An example is given by a pair of physicists who recently published their atomic scale study of the HIV-1 virus. Their simulation of a model containing more than 64 million atoms used force fields, making the simulation orders of magnitude cheaper than quantum mechanical calculations. Despite this enormous speedup, their simulation of 1.2 µs out of the life of an HIV-1 virus (actually it was only the outer skin of the virus, the inside was left empty) still took about 150 days on 3880 nodes of 16 cores on the Titan super computer of Oak Ridge National Laboratory (think about 25 512 years on your own computer).

In Flanders, scientist can make use of the TIER-1 facilities provided by the Flemish Super Computer (VSC). The first Tier-1 machine was installed and hosted at Ghent University. At the end of it’s life cycle the new Tier-1 machine (Breniac) was installed and is hosted at KULeuven. Although our Tier-1 supercomputer is rather modest compared to the Oak Ridge supercomputer (The HIV-1 calculation I mentioned earlier would require 1.5 years of full time calculations on the entire machine!) it allows Flemish scientists (including myself) to do things which are not possible on personal desktops or local clusters. I have been lucky, as all my applications for calculation time were successful (granting me between 1.5 and 2.5 million hours of CPU time every year). With the installment of the new supercomputer accounting of the requested resources has become fully integrated and automated. Several commands are available which provide accounting information, of which mam-balance is the most important one, as it tells how much credits are still available. However, if you are running many calculations you may want to know how many resources you are actually asking and using in real-time. For this reason, I wrote a small bash-script that collects the number of requested and used resources for running jobs:

13.         !$OMP END PARALLEL 14. 15. end subroutine OpenMPTest1 Notice in the example code above that outside the parallel section indicated with the directives$OMP PARALLEL and $OMP END PARALLEL, the program only sees a single thread, while inside the parallel section 8 threads will run (independent of the number of cores available). ### C. Simple parallelization The OpenMP frameworks consists of an set of directives which can be used to manage the parallelization of your code (cheat-sheets for OpenMP v3.0 and v4.0). I will not describe them in detail as there exists several very well written and full tutorials on the subject, we’ll just have a look at a quick and easy parallelization of a big for-loop. As said, OpenMP makes use of directives (or Pragma’s) which are placed as comments inside the code. As such they will not interfere with your code when it is compiled as a serial code (i.e. without the -fopenmp compiler flag). The directives are preceded by what is called a sentinel ($OMP ). In the above example code, we already saw a first directive: PARALLEL. Only inside blocks delimited by this directive, can your code be parallel.

 OpenMP second test
1. subroutine OMPTest2()
2.         use omp_lib;
3.
4.         integer :: IDT, NT,nrx,nry,nrz
5.         doubleprecision, allocatable :: A(:,:,:)
6.         doubleprecision :: RD(1:1000)
7.         doubleprecision :: startT, TTime, stt
8.
9.         call random_seed()
10.         call random_number(RD(1:1000))
11.         IDT=500 ! we will make a 500x500x500 matrix
12.         allocate(A(1:IDT,1:IDT,1:IDT))
13.
14.         write(*,'(A)') "Number of preferred threads:"
17.         startT=omp_get_wtime()
18.         !$OMP PARALLEL PRIVATE(stt) 19. stt=omp_get_wtime() 20. 21. !$OMP DO
22.         do nrz=1,IDT
23.            do nry=1,IDT
24.               do nrx=1,IDT
25.               A(nrx,nry,nrz)=RD(modulo(nrx+nry+nrz,1000)+1)
26.               end do
27.            end do
28.         end do
29.         !$OMP END DO 30. write(*,*) "time=",(omp_get_wtime()-stt)/omp_get_wtick()," ticks for thread ",omp_get_thread_num() 31. !$OMP END PARALLEL
32.         TTime=(omp_get_wtime()-startT)/omp_get_wtick()
33.         write(*,*)" CPU-resources:",Ttime," ticks."
34.
35.         deallocate(A)
36.     end subroutine RunTest2

The program above fills up a large 3D array with random values taken from a predetermined list. The user is asked to set the number of threads (lines 14-16), and the function omp_get_wtime() is used to obtain the number of seconds since epoch, while the function omp_get_wtick() gives the number of seconds between ticks. These functions can be used to get some timing data for each thread, but also for the entire program. For each thread, the starting time is stored in the variable stt. To protect this variable of being overwritten by each separate thread, this variable is declared as private to the thread (line 18: PRIVATE(stt) ). As a result, each thread will have it’s own private copy of the stt variable.

The DO directive on line 21, tells the compiler that the following loop needs to be parallelized. Putting the !$OMP DO pragma around the outer do-loop has the advantage that it minimizes the overhead produced by the parallelization (i.e. resources required to make local copies of variables, calculating the distribution of the workload over the different threads at the start of the loop, and combining the results at the end of the loop). As you can see, parallelizing a loop is rather simple. It takes only 4 additional comment lines (!$OMP PARALLEL , !$OMP DO, !$OMP END DO and !$OMP END PARALLEL) and some time figuring out which variables should be private for each thread, i.e. which are the variables that get updated during each cycle of a loop. Loop counters you can even ignore as these are by default considered private. In addition, the number of threads is set on another line giving us 5 new lines of code in total. It is of course possible to go much further, but this is the basis of what you generally need. Unfortunately, the presented example is not that computationally demanding, so it will be hard to see the full effect of the parallelization. Simply increasing the array size will not resolve this as you will quickly run out of memory. Only with more complex operations in the loop will you clearly see the parallelization. An example of a more complex piece of code is given below (it is part of the phonon-subroutine in HIVE):  Parallel example: Phonon spectra 1. !setup work space for lapack 2. N = this%DimDynMat 3. LWORK = 2*N - 1 4. call omp_set_num_threads(this%nthreads) 5. chunk=(this%nkz)/(this%nthreads*5) 6. chunk=max(chunk,1) 7. !$OMP PARALLEL PRIVATE(WORK, RWORK, DM, W, RPart,IO)
8.         allocate(DM(N,N))
9.         allocate( WORK(2*LWORK), RWORK(3*N-2), W(N) )
10.         !the write statement only needs to be done by a single thread, and the other threads do not need to wait for it
11.         !$OMP SINGLE 12. write(uni,'(A,I0,A)') " Loop over all ",this%nkpt," q-points." 13. !$OMP END SINGLE NOWAIT
14.         !we have to loop over all q-points
15.         !$OMP DO SCHEDULE(DYNAMIC,chunk) 16. do nrz=1,this%nkz 17. do nry=1,this%nky 18. do nrx=1,this%nkx 19. if (this%kpointListBZ1(nrx,nry,nrz)) then 20. !do nrk=1,this%nkpt 21. WORK = 0.0_R_double 22. RWORK = 0.0_R_double 23. DM(1:this%DimDynMat,1:this%DimDynMat)=this%dynmatFIpart(1:this%DimDynMat,1:this%DimDynMat) ! make a local copy 24. do nri=1,this%poscar%nrions 25. do nrj=1,this%poscar%nrions 26. Rpart=cmplx(0.0_R_double,0.0_R_double) 27. do ns=this%vilst(1,nri,nrj),this%vilst(2,nri,nrj) 28. Rpart=Rpart + exp(i*(dot_product(this%rvlst(1:3,ns),this%kpointList(:,nrx,nry,nrz)))) 29. end do 30. Rpart=Rpart/this%mult(nri,nrj) 31. DM(((nri-1)*3)+1:((nri-1)*3)+3,((nrj-1)*3)+1:((nrj-1)*3)+3) = & 32. & DM(((nri-1)*3)+1:((nri-1)*3)+3,((nrj-1)*3)+1:((nrj-1)*3)+3)*Rpart 33. end do 34. end do 35. call MatrixHermitianize(DM,IOS=IO) 36. call ZHEEV( 'V', 'U', N, DM, N, W, WORK, LWORK, RWORK, IO ) 37. this%FullPhonFreqList(:,nrx,nry,nrz)=sign(sqrt(abs(W)),W)*fac 38. end if 39. end do 40. end do 41. end do 42. !$OMP END DO
43.         !$OMP SINGLE 44. write(uni,'(A)') " Freeing lapack workspace." 45. !$OMP END SINGLE NOWAIT
46.         deallocate( WORK, RWORK,DM,W )
47.         !$OP END PARALLEL In the above code, a set of equations is solved using the LAPACK eigenvalue solver ZHEEV to obtain the energies of the phonon-modes in each point of the Brillouin zone. As the calculation of the eigenvalue spectrum for each point is independent of all other points, this is extremely well-suited for parallelization, so we can add !$OMP PARALLEL and !$OMP END PARALLEL on lines 7 and 47. Inside this parallel section there are several variables which are recycled for every grid point, so we will make them PRIVATE (cf. line 7, most of them are work-arrays for the ZHEEV subroutine). Lines 12 and 44 both contain a write-statement. Without further action, each thread will perform this write action, and we’ll end up with multiple copies of the same line (Although this will not break your code it will look very sloppy to any user of the code). To circumvent this problem we make use of the !$OMP SINGLE directive. This directive makes sure only 1 thread (the first to arrive) will perform the write action. Unfortunately, the SINGLE block will create an implicit barrier at which all other threads will wait. To prevent this from happening, the NOWAIT clause is added at the end of the block. In this specific case, the NOWAIT clause will have only very limited impact due to the location of the write-statements. But this need not always to be the case.

On line 15 the !\$OMP DO pragma indicates a loop will follow that should be parallelized. Again we choose for the outer loop, as to reduce the overhead due to the parallelization procedure. We also tell the compiler how the work should be distributed using the SCHEDULE(TYPE,CHUNK) clause. There are three types of scheduling:

1. STATIC: which is best suited for homogeneous workloads. The loop is split in equal pieces (size given by the optional parameter CHUNK, else equal pieces with size=total size/#threads)
2. DYNAMIC: which is better suited if the workload is not homogeneous.(in this case the central if-clause on line 19 complicates things). CHUNK can again be used to define the sizes of the workload blocks.
3. GUIDED: which is a bit like dynamic but with decreasing block-sizes.

From this real life example, it is again clear that OpenMP parallelization in fortran can be very simple.

### D. Speedup?

On my loyal sidekick (with hyper-threaded quad-core core i7) I was able to get following speedups for the phonon-code (the run was limited to performing only a phonon-DOS calculation):

Speedup of the entire phonon-subroutine due to parallelization of the main-phonon-DOS loop.

The above graph shows the speed-up results for the two different modes for calculating the phonon-DOS. The reduced mode (DM red), uses a spectrum reduced to that of a unit-cell, but needs a much denser sampling of the Brillouin zone (second approach), and is shown by the black line. The serial calculation in this specific case only took 96 seconds, and the maximum speedup obtained was about x1.84. The red and green curves give the speedup of the calculation mode which makes use of the super-cell spectrum (DM nored, i.e. much larger matrix to solve), and shows for increasing grid sizes a maximum speedup of x2.74 (serial time: 45 seconds) and x3.43 (serial time 395 seconds) respectively. The reason none of the setups reaches a speedup of 4 (or more) is twofold:

1. Amdahl’s law puts an upper limit to the global speedup of a calculation by taking into account that only part of the code is parallelized (e.g. write sections to a single file can not be parallelized.)
2. There needs to be sufficient blocks of work for all threads (indicated by nkz in the plot)

In case of the DM nored calculations, the parallelized loop clearly takes the biggest part of the calculation-time, while for the DM red calculation, also the section generating the q-point grid takes a large fraction of the calculation time limiting the effect of parallelization. An improvement here would be to also parallelize  the subroutine generating the grid, but that will be for future work. For now, the expensive DM nored calculations show an acceptable speedup.

## Folding Phonons

About a year ago, I discussed the possibility of calculating phonons (the collective vibration of atoms) in the entire Brillouin zone for Metal-Organic Frameworks. Now, one year later, I return to this topic, but this time the subject matter is diamond. In contrast to Metal-Organic Frameworks, the unit-cell of diamond is very small (only 2 atoms). Because a phonon spectrum is calculated through the gradients of forces felt by one atom due to all other atoms, it is clear that within one diamond unit-cell these forces will not be converged. As such, a supercell will be needed to make sure the contribution, due to the most distant atoms, to the experienced forces, are negligible.

Using such a supercell has the unfortunate drawback that the dynamical matrix (which is $3N \times 3N$, for N atoms) explodes in size, and, more importantly, that the number of eigenvalues, or phonon-frequencies also increases ($3N$) where we only want to have 6 frequencies ( $3 \times 2$ atoms) for diamond. For an $M \times M \times M$ supercell we end up with $24M^3 -6$  additional phonon bands which are the result of band-folding. Or put differently, $24M^3 -6$ phonon bands coming from the other unit-cells in the supercell. This is not a problem when calculating the phonon density of states. It is, however, a problem when one is interested in the phonon band structure.

The phonon spectrum at a specific q-point in the first Brillouin zone is given by the square root of the eigenvalues of the dynamical matrix of the system. For simplicity, we first assume a finite system of n atoms (a molecule). In that case, the first Brillouin zone is reduced to a single point q=(0,0,0) and the dynamical matrix looks more or less like the hessian:

With $\varphi (N_a,N_b) = [\varphi_{i,j}(N_a , N_b)]$ $3 \times 3$ matrices $\varphi_{i,j}(N_a,N_b)=\frac{\partial^2\varphi}{\partial x_i(N_a) \partial x_j(N_b)} = - \frac{\partial F_i (N_a)}{\partial x_j (N_b)}$  with i, j = x, y, z. Or in words, $\varphi_{i,j}(N_a , N_b)$ represents the derivative of the force felt by atom $N_a$ due to the displacement of atom $N_b$. Due to Newton’s second law, the dynamical matrix is expected to be symmetric.

When the system under study is no longer a molecule or a finite cluster, but an infinite solid, things get a bit more complicated. For such a solid, we only consider the symmetry in-equivalent atoms (in practice this is often a unit-cell). Because the first Brillouin zone is no longer a single point, one needs to sample multiple different points to get the phonon density-of-states. The role of the q-point is introduced in the dynamical matrix through a factor $e^{iq \cdot (r_{N_a} - r_{N_b}) }$, creating a dynamical matrix for a single unit-cell containing n atoms:

Because a real solid contains more than a single unit-cell, one should also take into account the interactions of the atoms of one unit-cell with those of all other unit-cells in the system, and as such the dynamical matrix becomes a sum of matrices like the one above:

Where the sum runs over all unit-cells in the system, and Ni indicates an atom in a specific reference unit-cell, and MRi  an atom in the Rth unit-cells, for which we give index 1 to the reference unit-cell. As the forces decay with the distance between the atoms, the infinite sum can be truncated. For a Metal-Organic Framework a unit-cell will quite often suffice. For diamond, however, a larger cell is needed.

An interesting aspect to the dynamical matrix above is that all matrix-elements for a sum over n unit-cells are also present in a single dynamical matrix for a supercell containing these n unit-cells. It becomes even more interesting if one notices that due to translational symmetry one does not need to calculate all elements of the entire supercell dynamical matrix to construct the full supercell dynamical matrix.

Assume a 2D 2×2 supercell with only a single atom present, which we represent as in the figure on the right. A single periodic copy of the supercell is added in each direction. The dynamical matrix for the supercell can now be constructed as follows: Calculate the elements of the first column (i.e. the gradient of the force felt by the atom in the reference unit-cell, in black, due to the atoms in each of the unit-cells in the supercell). Due to Newton’s third law (action = reaction), this first column and row will have the same elements (middle panel).

Translational symmetry on the other hand will allow us to determine all other elements. The most simple are the diagonal elements, which represent the self-interaction (so all are black squares). The other you can just as easily determine by looking at the schematic representation of the supercell under periodic boundary conditions. For example, to find the derivative of the force on the second cell (=second column, green square in supercell) due to the third cell (third row, blue square in supercell), we look at the square in the same relative position of the blue square to the green square, when starting from the black square: which is the red square (If you read this a couple of times it will start to make sense). Like this, the dynamical matrix of the entire supercell can be constructed.

This final supercell dynamical matrix can, with the same ease, be folded back into the sum of unit-cell dynamical matrices (it becomes an extended lookup-table). The resulting unit-cell dynamical matrix can then be used to create a band structure, which in my case was nicely converged for a 4x4x4 supercell. The bandstructure along high symmetry lines is shown below, but remember that these are actually 3D surfaces. A nice video of the evolution of the first acoustic band (i.e. lowest band) as function of its energy can be found here.

The phonon density of states can also be obtained in two ways, which should, in contrast to the band structure, give the exact same result: (for an $M \times M \times M$ supercell with n atoms per unit-cell)

1. Generate the density of states for the supercell and corresponding Brillouin zone. This has the advantage that the smaller Brillouin zone can be sampled with fewer q-points, as each q-point acts as M3 q-points in a unit-cell-approach. The drawback here is the fact that for each q-point a (3nM3)x(3nM3) dynamical matrix needs to be solved. This solution scales approximately as O(N3) ~ (3nM3)3 =(3n)3M9. Using linear algebra packages such as LAPACK, this may be done slightly more efficient (but you will not get O(N2) for example).
2. Generate the density of states for the unit-cell and corresponding Brillouin zone. In this approach, the dynamical matrix to solve is more complex to construct (due to the sum which needs to be taken) but much smaller: 3nx3n. However to get the same q-point density, you will need to calculate M3 times as many q-points as for the supercell.

In the end, the choice will be based on whether you are limited by the accessible memory (when running a 32-bit application, the number of q-point will be detrimental) or CPU-time (solving the dynamical matrix quickly becomes very expensive).