Large systems: the Fullerene molecule

In this tutorial, by solving the Fullerene molecule, we show you how to do simulations of systems with a large number of atoms with Octopus. Most of the default parameters work fine for large systems, but there are a few things you might need to adjust.

The eigensolver

The default eigensolver of Octopus is simple and reliable, but for large systems (> 50 orbitals) it can be quite slow. For these systems it is better to use the RMMDIIS eigensolver. This is a very efficient diagonalizer, but it can have convergence problems. In Octopus we follow the algorithm implementation of Kresse and Furthmüller 1 that provides a very robust algorithm if the following considerations are taken:

Let’s see an example: this is an input file, inp , to calculate the ground state of the fullerene molecule with the RMMDIIS eigensolver:


 XYZCoordinates = "C60.xyz"
 UnitsOutput = eV_Angstrom
 CalculationMode = gs
 EigenSolver = rmmdiis
 ExtraStates = 20

In this case, to use RMMDIIS the only thing we need to do is to change the Eigensolver variable and to include 20 ExtraStates (by default Octopus does not include any unoccupied states).

This is the content of the coordinates file, C60.xyz :


60

 C    -1.415441     3.015027    -1.167954
 C     2.583693     2.292084    -0.721147
 C     0.721155    -2.583700     2.292090
 C     3.015024    -1.167941     1.415429
 C    -3.015019    -1.167950     1.415434
 C    -0.721134    -2.583697     2.292085
 C    -2.583698     2.292079    -0.721148
 C    -3.461397    -0.000006     0.694328
 C     1.415423     3.015023    -1.167951
 C     3.461399     0.000011     0.694319
 C    -1.415430    -3.015023    -1.167948
 C    -2.292086    -0.721156    -2.583690
 C    -0.000003     0.694321    -3.461398
 C     2.292087    -0.721146    -2.583706
 C     1.415436    -3.015019    -1.167948
 C    -2.292087     0.721141    -2.583692
 C     1.167941     1.415428    -3.015025
 C     3.015020    -1.167945    -1.415437
 C     0.694330    -3.461397     0.000006
 C    -2.583701    -2.292092    -0.721145
 C     0.721144     2.583697     2.292082
 C     1.167948     1.415436     3.015018
 C    -0.000002     0.694327     3.461397
 C    -1.167951     1.415435     3.015031
 C    -0.721152     2.583693     2.292084
 C    -0.694329     3.461397    -0.000006
 C     2.583701     2.292093     0.721144
 C     2.292088    -0.721140     2.583692
 C    -1.167941    -1.415428     3.015025
 C    -3.015020     1.167945     1.415437
 C     1.167955    -1.415429     3.015019
 C    -2.292077    -0.721143     2.583701
 C    -2.583697     2.292080     0.721148
 C     0.694318     3.461398    -0.000005
 C     3.015030     1.167958     1.415430
 C    -0.694317    -3.461398     0.000005
 C    -3.015030    -1.167958    -1.415430
 C    -1.167955     1.415429    -3.015019
 C     2.292078     0.721143    -2.583701
 C     2.583697    -2.292079    -0.721148
 C     1.167951    -1.415434    -3.015031
 C     0.721152    -2.583693    -2.292084
 C    -0.721144    -2.583696    -2.292082
 C    -1.167948    -1.415436    -3.015019
 C     0.000002    -0.694327    -3.461397
 C     0.721134     2.583697    -2.292085
 C     3.461397     0.000007    -0.694329
 C     1.415440    -3.015027     1.167954
 C    -2.583694    -2.292084     0.721147
 C    -3.015025     1.167942    -1.415429
 C    -2.292087     0.721146     2.583706
 C    -1.415436     3.015019     1.167948
 C     1.415430     3.015023     1.167947
 C     2.292086     0.721156     2.583690
 C     0.000003    -0.694321     3.461398
 C    -1.415423    -3.015023     1.167952
 C    -3.461399    -0.000011    -0.694320
 C    -0.721155     2.583700    -2.292090
 C     3.015019     1.167950    -1.415434
 C     2.583698    -2.292079     0.721148

Create both files and run Octopus. It should take a few minutes (depending on the speed of your machine). There are some peculiarities of the RMMDIIS eigensolver that you might notice:

Self-consistency convergence

Once your calculation is finished, you might notice that the final residue in the eigenvectors is rather large (of the order of $10^{-4}$). This is for two reasons: first, it is possible to converge the density without properly converging the eigenvectors sometimes. Try a tighter tolerance and rerun (with restart), adding


 ConvRelDens = 1e-6

to your inp file. Now your eigenvectors should be better converged.

Can you converge all the eigenvectors? How many states are actually properly converged? Suppose you wanted to calculate correct unoccupied states. How would you do it and be sure they were converged?

The second reason for poor convergence is that we are using a large (default) spacing – check what it is in the output. The problem is just not as well-conditioned with a large spacing. So, you can try using a smaller spacing such as the one you found in the methane tutorial.