# Difference between revisions of "Tutorial:Total energy convergence"

(→Input) |
|||

Line 89: | Line 89: | ||

{{tutorial_foot|next=Tutorial:Centering geometry|prev=Tutorial:Nitrogen atom}} | {{tutorial_foot|next=Tutorial:Centering geometry|prev=Tutorial:Nitrogen atom}} | ||

+ | |||

+ | [[Category:Basic|M]] |

## Revision as of 15:18, 17 July 2018

In this example, for the methane molecule CH_{4}, we will play a little bit with the spacing and the dimensions of the mesh.

## Input

In the input file we define a variable CH that represents the bond length between the carbon and the hydrogen atoms. From our basic chemistry class we know that methane has a tetrahedral structure. If we put the carbon atom at the origin, the hydrogen atoms have the coordinates given in the following input file.

`CalculationMode`

= gs`UnitsOutput`

= eV_Angstrom`Radius`

= 3.5*angstrom`Spacing`

= 0.22*angstrom CH = 1.2*angstrom %`Coordinates`

"C" | 0 | 0 | 0 "H" | CH/sqrt(3) | CH/sqrt(3) | CH/sqrt(3) "H" | -CH/sqrt(3) |-CH/sqrt(3) | CH/sqrt(3) "H" | CH/sqrt(3) |-CH/sqrt(3) | -CH/sqrt(3) "H" | -CH/sqrt(3) | CH/sqrt(3) | -CH/sqrt(3) %

We start with a bond length of 1.2 Å but that's not so important at the moment, since we can optimize it later (for more information see the tutorial on geometry optimization). (You should not use 5 Å or so, but something slightly bigger than 1 Å is fine.)

Some notes concerning the input file:

- We do not tell Octopus explicitly which
`BoxShape`

to use. The default is a union of spheres centered around each atom. This turns out to be the most economical choice in almost all cases. - We also do not specify the %
`Species`

block. In this way, Octopus will use default pseudopotentials for both Carbon and Hydrogen. This should be OK in many cases, but, as a rule of thumb, you should do careful testing before using any pseudopotential for serious calculations.

## Convergence

If you use the given input file you should find the following values in the resulting `static/info` file.

Eigenvalues [eV] #st Spin Eigenvalue Occupation 1 -- -15.988934 2.000000 2 -- -9.064039 2.000000 3 -- -9.064039 2.000000 4 -- -9.064039 2.000000 Energy [eV]: Total = -219.01537542 Free = -219.01537542 ----------- Ion-ion = 236.08498119 Eigenvalues = -86.36210292 Hartree = 393.44961913 Int[n*v_xc] = -105.38115372 Exchange = -69.66918783 Correlation = -11.00060043 vanderWaals = 0.00000000 Delta XC = 0.00000000 Entropy = 0.00000000 -TS = -0.00000000 Kinetic = 163.96741296 External = -931.84569058 Non-local = -49.74424172

Now the question is whether these values are converged or not. This will depend on the `Spacing`

and the `Radius`

. The only way to answer this question is to try other values for these variables. We will start with the spacing. Since we are interested in the total energy of the system, we will look at the changes of its value. We will keep all entries in the input file fixed except for the spacing that we will make smaller by 0.02 Å all the way down to 0.1 Å. So you have to run Octopus several times (you can use the script from the Nitrogen tutorial) to get the following results.

#Spacing Total Energy 0.22 -219.01537542 0.20 -218.58289833 0.18 -218.20446230 0.16 -218.14597804 0.14 -218.14033248 0.12 -218.13672950 0.10 -218.12988174

If you give these numbers to gnuplot (or other software to plot) you will get a curve like the one shown on the right.

As you can see from this picture, the total energy is converged to within 0.1 eV for a spacing of 0.18 Å. So we will use this spacing for the next calculations, and play with the `Radius`

of the box. We will change the radius in steps of 0.5 Å. Doing this we get

#Radius Total Energy 2.5 -217.99419597 3.0 -218.16990643 3.5 -218.20446230 4.0 -218.21157996 4.5 -218.21300375 5.0 -218.21356900

If we again ask for a convergence up to 0.1 eV we should use a radius of 3.5 Å. How does this compare to the size of the molecule? Can you explain why the energy increases (becomes less negative) when one decreases the size of the box?