The Gaussian Plane Waves method (GPW) solves the DFT Kohn-Sham equations efficiently. It uses gaussians as basisset, and planewaves as auxiliary basis. This is similar at the Resolution of Identity (RI) methods but with a different basisset.
In GPW the whole density is transferred to plane waves, and one has the density $n(r)=\sum_{i j} P_{i j} \phi_i (r) \phi_j(r)$ in the gaussian basis set and the density $\tilde n$ taken on a uniform grid.
Then $\tilde n$ is used to calculate the hartree potential $V_H$ (via an FFT based poisson solver) and the exchange and correlation potential $V_\text{xc}$. These potential are then transferred back to gaussian basis set by integrating them with $\phi_i(r)\phi_j(r)$. To make the collocation and integration perfectly consistent with each other one can set FORCE_EVAL%DFT%QS%MAP_CONSISTENT this normally adds only a very small overhead to the calculation and is
$n$ and $\tilde n$ are not equal, and this introduces an error in the calculation. $\tilde n$ converges toward $n$ when the cutoff (that controls the grid spacing) goes to infinity (and gridspacing to 0). Which cutoff is sufficient to represent a density depends on how sharp is the gaussian basis set (or that of the potential, but it is always broader).
For historical reasons the density of the grid is given as the energy (in Ry) of the highest reciprocal vector that can be represented on the grid. This can be roughly given as $0.5(\pi/dr)^2$ where $dr$ is the gridspacing in Bohr. The characteristic length of a gaussian with exponent A is given by $1/\sqrt{a}$ (up to a factor 2 depending on the convention used). This means that the cutoff to represent in the same way a gaussian depends linearly on the exponent. Thus one can get a first guess for an acceptable guess can be take from the knowledge that for water with $\alpha_H=47.7$ a good cutoff for doing MD is 280 Ry.
It turns out that if one wants to put the whole density on the grid, the core electrons of even the simplest atoms cannot be represented, thus one has to remove to core electrons and use pseudopotentials for the atoms. In CP2K we use the Goedecker-Teter-Hutter pseudopotentials, these are harder than other pseudopotentials, but also more accurate.
$\tilde n$ is optimized for the electrostatic part, but is used also to calculate the exchange and correlation potential. Because of this, and because the GTH pseudopotential goes almost to 0 close to the core of the atom, the xc potential, especially for gradient corrected functionals, converges badly. Instead of using very high cutoffs one can perform a smoothing of the density, and calculate the derivatives on the grid with other methods than the G-space based derivatives.
For MD of water using a cutoff of 280 Ry XC_SMOOTH_RHO NN10 and XC_DERIV SPLINE2_SMOOTH (in the FORCE_EVAL%DFT%XC%XC_GRID section) give good results, please note that these options renormalize the total energy, and the amount of renormalization is dependent on the cutoff. Thus energies with different cutoffs cannot be easily compared, only interaction energies or forces can be calculated.
Methods that do not redefine the total energy are XC_SMOOTH_RHO NONE and XC_DERIV equal to either PW, SPLINE3 or SPLINE2. These are listed from the one that assumes more regularity (PW the the one that assumes less regularity SPLINE2. Normally SPLINE2 is a good choice, but for high cutoffs (600 Ry for water) SPLINE3 is better. The default (PW) is not bad, but generally inferior to the others.
The QS part of Cp2k uses basis sets that are a linear combination gaussian functions. As the GPW method uses pseudopotentials one cannot use basis set of other programs, at least the core part should be optimized for cp2k. The polarization and augmentation functions can be taken from dunning type basis sets.
Cp2k basis normally are build with an increasing accuracy level starting from SZV (single Z valence, normally only for really crude results), DZVP (double Z valence and one polarization function, already suitable for MD, and give reasonable results), TZV2P (triple Z valence and two polarization functions, fairly good results), QZV3P (quadruple Z valence and three polarization functions, good results). Note that the discussion about the quality of the result is very indicative, and valid for condensed phase MD, for gas phase to reduce the BSSE an augmented basis (aug-) with diffuse basis functions is needed, if you calculate properties that depend on the polarization probably also (aug-) and polarization functions will be important. In any case you should perform a convergence check for your properties with respect to the basis set.
A good set of basis set are available in the BASIS_SET and GTH_BASIS_SETS files. Other can be created with a program, and you can ask on the list or try to optimize them with the program that is part of cp2k.
The basis set to use for an element is set in the FORCE_EVAL%SUBSYS%KIND section with its name with the BASIS_SET keyword.
The GTH pseudopotential that one has to use has a large library of elements available in the POTENTIAL file. The pseudopotential has to be optimized for the exchange correlation functional that one is using. Normally changing from one functional to the other is not too difficult with the optimization program that is part of cp2k. If your element-functional combination is missing ask on the forum.
The pseudopotential to use for an element is set in the FORCE_EVAL%SUBSYS%KIND section with its name.
In DFT the choice exchange correlation functional is an important issue, because unfortunately results can depend on it. Here we don't want to discuss how to select it, normally it is a good idea to use the same one as what is used in your field, so that you can more easily compare your results with the literature. An important choice though is whether to use exact exchange or an hybrid functional or not. Hybrid functionals normally improve the description of barriers and radicals, but are much more costly. In general one should start with a GGA and meta GGA and switch to hybrid functionals to check the results or if needed.
In cp2k the xc functional is selected in the FORCE_EVAL%DFT%TDDFPT%XC section. Common functional combinations can be selected directly as section parameters of the XC_FUNCTIONAL for example with the POTENTIAL keyword.</p>
&XC_FUNCTIONAL BLYP &END XC_FUNCTIONAL
but more complex combination have to be chosen by directly setting the underlying subsections.
Apart from cutoff and basis set with the GPW method there are two other important switches that control the convergence of the method. One is FORCE_EVAL%DFT%QS%EPS_DEFAULT that controls the convergence of integral calculation, overlap cutoff, neighboring lists… and sets a slew of other EPS_* keywords so that the KS matrix is calculated with the requested precision. The other is EPS_SCF which controls the convergence of the scf procedure. In general one should have an error $e_f$ on the forces one should have $\sqrt{\text{eps}_\text{scf}} < e_f $ and $\sqrt{\text{eps}_\text{default}} < \text{eps}_\text{scf}$, and in general around 1.e-3 is the error that you need to have reasonable MD.
To perform the SCF and find the groundstate there are two main methods: the traditional diagonalization method accelerated by the DIIS procedure, and the orbital transform method.
FORCE_EVAL%DFT%SCF%SCF_GUESS RESTART is a parameter that controls how the initial density is built, normally one uses either ATOMIC or RESTART. During an MD this setting normally influences only the initial startup.
Traditional diagonalization is like most other DFT codes, it diagonalizes the KS matrix to get a new density from the n lowest eigenvectors.
The density for the next scf cycle is built at the beginning just mixing the new and the old density with the FORCE_EVAL%DFT%SCF%MIXING factor. If the new density is close enough to the old density FORCE_EVAL%DFT%SCF%EPS_DIIS then the DIIS procedure is activated.
Diagonalization works well, but for difficult or big systems the OT method is better (and in general there is no reason not to use it as default). The OT method directly minimizes the the electronic energy with respect to the wavefunctions. It uses a clever parametrization of the wavefunctions so that the orthogonality constraint becomes a linear constraint.
To activate OT adding the section FORCE_EVAL%DFT%SCF%OT is enough.
The advantage of the OT method, are that being a direct method it always converges to something, that it needs less memory than the diagonalization (important for big systems), and each SCF iteration is faster (but it might need a little more iterations). Anyway normally it is faster than diagonalization.
The speed of convergence of the OT method the preconditioner choice (FORCE_EVAL%DFT%SCF%PRECONDITIONER) is important. FULL_KINETIC is a good choice for very large systems, but for smaller or difficult systems FULL_ALL is normally better. The quality of the preconditioner is dependent on how close one is to the solution. Especially with expensive preconditioners like FULL_ALL being closer might improve much the preconditioner and the convergence. For this reason it might be worthwhile to rebuild the preconditioner after some iterations of the SCF. This might be reached using setting FORCE_EVAL%DFT%SCF%MAX_SCF to a small value, and activating the outer scf loop addint the FORCE_EVAL%DFT%SCF%OUTER_SCF section, and setting its EPS_SCF equal to the one of the scf loop, and chooing a outer MAX_SCF big enough (the total number of iteration is the product of internal and external MAX_SCF.
Also the minimizer of the OT method can be changed with FORCE_EVAL%DFT%SCF%OT%MINIMIZER, the default is CG that uses a very robust conjugated gradient method. Less robust but often faster is the DIIS method.