Difference between revisions of "Tutorial:RDMFT"
m (Micael moved page Category:RDMFT to Tutorial:RDMFT without leaving a redirect) 

(No difference)

Revision as of 12:46, 19 December 2019
In this tutorial, you will learn how to do a Reduced Density Matrix Functional Theory (RDMFT) calculation with Octopus.
In contrast to density functional theory or HartreeFock theory, here we do not try to find one optimal slaterdeterminant (singlereference) to describe the groundstate of a given manybody system, but instead approximate the onebody reduced density matrix (1RDM) of the system. Thus, the outcome of an RDMFT minimization is a set of eigenvectors, so called natural orbitals (NOs,) and eigenvalues, socalled natural occupation numbers (NONs), of the 1RDM. One importance aspect of this is that we need to use more' orbitals than the number of electrons. This additional freedom allows to include static correlation in the description of the systems and hence, we can describe settings or processes that are notoriously difficult with any singlereference method. One such process is the dissociation of a molecule, which thus becomes our example for this tutorial: You will learn how to do a properly converged dissociation study of H_2. A good recap of RDMFT can be found in ...
Basis Steps of an RDMFT Calculation
The first thing to notice for RDMFT calculations is that, we will explicitly make use of a basis set, which is in contrast to the minimization over the full realspace grid that is performed normally with Octopus. So we always need to do a preliminary calculation to generate our basis set. So let us start with the hydrogen molecule in equilibrium position. Create the following inp file:
CalculationMode
= gsTheoryLevel
= independent_particlesDimensions
= 3Radius
= 8Spacing
= 0.15 # distance between H atoms d = 1.4172 # (equilibrium bond length H2) %Coordinates
'H'  0  0  d/2 'H'  0  0  d/2 %ExtraStates
= 14
You should be already familiar with all these variables, so there is no need to explain them again, but we want to stress the two important points here:
 We chose
TheoryLevel
= independent_particles, which you should always use as a basis in an RDMFT calculation (see actual octopus paper?)  We set the
ExtraStates
variable, which controls the size of the basis set that will be used in the RDMFT calculation. Thus in RDMFT, we have with the basis set a new numerical parameter besides theRadius
and theSpacing
that needs to be converged for a successful calculation. The size of the basis set M is just the number of the orbitals for the ground state (number of electrons divided by 2) plus the number of extra states.  All the numerical parameters depend on each other: If we want to include many ExtraStates to have a sufficiently large basis, we will need also a larger Radius and especially a smaller Spacing at a certain point. The reason is that the additional states will have function values bigger than zero in a larger region than the bound states. Additionally all states are orthogonal to each other, and thus the grid needs to resolve more and more nodes.
With this basis, we can now do our first rdmft calculation. For that, you just need to change the TheoryLevel
to 'rdmft' and add the part ExperimentalFeatures
= yes. In the standard out there will be some new information:
 Calculating Coulomb and exchange matrix elements in basis this may take a while
 Occupation numbers: they are not zero or one/two here but instead fractional!
 anything else important?
Basis Set Convergence
To find good value for the Radius and Spacing, one should perform convergence sets with a larger number of ExtraStates, in our case say 14 (thus M=1+19=20) and make sure that they are converged. Having done this, we can go on and converge the basis set. So let us perform a series.
For that, we first create large enough basis by repeating the above calculation with say ExtraStates = 29. We copy the restart folder in a new folder that is called 'basis' and execute the following script (if you rename anything you need to change the respective part in the script):
#!/bin/bash series=ESseries outfile="$series.log" echo "#ES Energy 1. NON 2.NON" > $outfile list="4 9 14 19 24 29" export OCT_PARSE_ENV=1 for param in $list do folder="$series$param" mkdir $folder out_tmp=outRDMFT$param cd $folder # here we specify the basis folder cp r ../basis/restart . # create inp file { cat <<EOF calculationMode = gs TheoryLevel = rdmft ExperimentalFeatures = yes ExtraStates = $param Dimensions = 3 Radius = 8 Spacing = 0.15 # distance between H atoms d = 1.4172 # (equilibrium) %Coordinates 'H'  0  0  d/2 'H'  0  0  d/2 % Output = density OutputFormat = axis_x EOF } > inp # put the octopus dir here [octopusdir]/bin/octopus > $out_tmp energy=`grep a "Total energy" $out_tmp  tail 1  awk '{print $3}'` seigen=`grep a " 1 " static/info  awk '{print $2}'` peigen=`grep a " 2 " static/info  awk '{print $2}'` echo $param $energy $seigen $peigen >> ../$outfile cd .. done
If everything works out fine, you should find the following values in the 'ESseries.log' file:
#ES Energy 1. NON 2.NON 4 1.1476150018E+00 1.935750757008 0.032396191164 9 1.1498006205E+00 1.932619193929 0.032218954381 14 1.1609241676E+00 1.935215440985 0.032526426664 19 1.1610006378E+00 1.934832116929 0.032587169713 24 1.1622104536E+00 1.932699204653 0.032997371081 29 1.1630839347E+00 1.932088486112 0.032929702131
So we see that even with M=30, the calculation is not entirely converged! Thus for a correctly converged result, one would need to further increase the basis and probably resort to a cluster calculation. For the purposes of this tutorial, we will use not entirely converged parameters.
H2 Dissociation
Now, we can go on with our original goal: the H2 dissociation curve.