This document is thought as way to get you started with cp2k. It does not focus on describing how to calculate what you are interested in, your properties, but it should help you to perform meaningful simulations, avoiding stupid technical errors and without too much overhead. The basic structure of this comes from an informal tutorial on cp2k that I gave in Prag. If you have suggestions or improvements for it you can contact me at fawzi@gmx.ch. If you have further questions about cp2k, then go to http://cp2k.berlios.de/ and ask in the forum, normally the developers will be glad to help you out, and maybe if you want, also to start a collaboration with you.
One of the easiest ways to get cp2k if one has cvs installed (and it is a good idea) is
touch ~/.cvspass
You can add
cvs -d :pserver:anonymous@cvs.cp2k.berlios.de:/cvsroot/cp2k login
cvs -d :pserver:anonymous@cvs.cp2k.berlios.de:/cvsroot/cp2k co cp2k
-D YYYY-MM-DDTHH:MM to specify the date.
To compile it you have to adapt the architecture file in cp2k/arch. Executing cp2k/tools/get_arch_name you can see your architecture, changing the environment variable FORT_C_NAME you can change the compiler on some architectures.
As make you should use the gnu make, and the architecture can be changed also with make ARCH=my_arch_name.
The kind of build (sdbg, sopt, pdbg, popt) can be changed with VERSION=build_type (actually using the build type as target also works).
Targets can be: build, lib, clean, realclean, distclean.
CP2K has an extensive series of test inputs that are executed automatically and it is a good idea to test your executable, this way you can be quite confident that the executable you have can be trusted. In fact CP2K uses a lot of complier features that can expose bugs of the compiler or incorrect built libraries.
To setup the regtests you have to
cp2k directory and of cp2k/tools/do_regtest into it for example withmkdir tst
cd tst
cvs -d :pserver:anonymous@cvs.cp2k.berlios.de:/cvsroot/cp2k co cp2k
cp cp2k/tools/do_regtest .
cp2k and do_regtest), decide the version you want to check and set dir_triplet (normally to what is returned by cp2k/tools/get_arch_name).
you probably also
export FORT_C_NAME=g95
dir_base=/home/regtest/autotest/g95
cp2k_version=sopt
dir_triplet=Linux-x86-64-g95
it is commented.
Sometime for it to work you need to define datum like
date -u '+%Y-%m-%dT%H:%M:%S' to avoid empty spaces in it.
Documentation of cp2k can be found on the website http://cp2k.berlios.de/ and can also be generated by the cp2k executable itself using the flag --html or --xml (followed by xml2htm if you want an html documentation).
The documentation generated by the executable is always consistent with its version, whereas the one on the website refers to the latest version.
KEYWORD FALSE).OT section has a logical default parameter that controls the activation of the OT method. So,
&OT
...
&END OT&OT FALSE
...
&END OTThe basic idea of the input is to define the options for the method that you want to use to calculate the forces/energy in the FORCE_EVAL section.
If you use the ab-initio DFT part then the DFT and QS subsections will probabily be of interest for you.
The method to use is choosen with FORCE_EVAL%METHOD keyword, i.e. the keyword METHOD in the FORCE_EVAL section. So you can easily switch from one method to the another if you gave meaningful parameters to both.
Then parameters for the kind of calculation (run type) to do using that method (energy,force,molecular dynamics,langevin,...) are defined in the MOTION section
Finally in the GLOBAL section you can choose the run type to perform with the keyword RUN_TYPE, the name of the project (PROJECT), used to generate most of the output files, the PRINT_LEVEL keyword that controls the general level of output, and other general flags. The GLOBAL%PRINT%PHYSCON keyword enables a printout of the physical constants used by cp2k.
To control the output of information or properties often printkeys are used. A printkey is a section that has a default value that says at which print level to start printing the corresponding values, and keywords to control their output. Normally printkeys are kept together in a PRINT section.
To understand how to control the output it is necessary to first understand the iteration level. During the calculations cp2k keeps a global iteration level which is a set of integers that keep track which iteration one is doing at that moment.
The exact meaning of each iteration level depends on the run type, for example doing an MD with the QS GPW method an iteration level of 1_4_5 means that the program is doing the 5th scf step of the 4th md step, whereas 1_9 means the we are at the 9th md setp (outside the scf loop). The first 1 is a root iteration level that stays always 1 during a normal md run
| section parameter | Level starting at which this property is printed, valid values are
ON, OFF, SILENT, LOW, MEDIUM, HIGH, DEBUG
No section means the default value, the presence of the prind sections implies a value of ON (always printed), if one wants to selectively switch off the printkey it is possible the give it the explicit value OFF, i.e. &SECTION OFF
|
| EACH | How often this proprety is printed, this is matched with the actual iteration level from the right replacing non present levels with 1. For example "2 1" for an scf property during an scf means to write each scf step each 2 md step. How to handle the last iteration is treated separately in ADD_LAST (this mean that EACH 0 might print the last iteration) |
| ADD_LAST | If the last iteration should be added, and if it should be marked symbolically (with l) or with the iteration number. Not every process is able to identify the last iteration early enough to be able to output. Valid keywords: NO, NUMERIC, SYMBOLIC |
| COMMON_ITERATION_LEVELS | How many iterations levels should be written in the same file (no extra information about the actual iteration level is written to the file). If it is 0 then all iterations go to different files, if it is 1 then (if the current property is an SCF property and one is doing an MD then each MD step goes to a separate file (but different scf steps are in the same file) |
| FILENAME | controls part of the filename for output.
|
The files to restart a trajectory in cp2k are just input files with all the defaults made explicit. So it should not be so difficult to understand them.
Often one whats to continue a trajectory with slightly different parameters, for example of an equilibration run followed by a NVE run, and having to edit the restart file is not so practical, and cannot be prepared before the previous run has finished. To make this easier one can prepare the input that takes part of its values from another input. In particular it might take positions and velocities from an external input.
This is done with the section EXT_RESTART giving the name of the external input with the EXTERNAL_FILE_NAME keyword. In thsi section there are keywords that controls exactly what is taken from the external file, for example if one wants to restart just the positions and not the velocities he should set RESTART_VEL FALSE.
In DFT calculation there is also another restart, the one of the wavefunctions. This restart is stored in a binary file that is architecture specific (see cp2k/tools/RestartTools if you need to transfer it). The name and location of this file are controlled by FORCE_EVAL%DFT%RESTART_FILE_NAME.
Now that you understand the main structure of an input you should be able to better understand the sample files that are present in the cp2k/test directory. Please note that only the files used for regtest are guaranteed to stay up to date, but the regtest are optimized to have a small running time, not to have meaningful results, so be careful in using them blindly.
The Gaussian Plane Waves method [Lip97] to efficiently solve the DFT Kohn-Sham equations. 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)=sumi jPi j phii(r)phij(r) in the gaussian basis set and the density ñ taken on a uniform grid.
Then ñ is used to caluclate the hartree potential VH (via an FFT based poisson solver) and the exchange and correlation potential Vxc. These potential are then transferred back to gaussian basis set by integrating them with phii(r)phij(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 ñ are not equal, and this introduces an error in the calculation. ñ 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 hystorical reasons the density of the grid is given as the energy (in Ry) of the highest reciprocal vector that can be represented on the 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 alphaH=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 Godeker-Theta-Hutter pseudopotentials [GTH96], these are harder than other pseudopotentials, but also more accurate.
ñ is optimized for the electrostatic part, but is used also to calculate the echange and correlation potential. Because of this, and because the GTH pseudopotential goes almost to 0 close to the core of the atom, the xc potentisl, 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%SC%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. Thes 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.
&XC_FUNCTIONAL BLIP
&END XC_FUNCTIONALApart from cutoff and basis set with the GPW method there are two other importan 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 ef on the forces one should have sqrt(eps_scf) < ef and sqrt(eps_default) < eps_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 disgonalization method accellerated 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 [VH03] 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 parametrizations 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
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 withFORCE_EVAL%DFT%SCF%OT%MINIMIZER, the default is CG that uses a very robust conjugated gradient method. Less roubust but often faste is the DIIS method.
Molecular dynamics a good method to perform themodynamical averages, and to look at dynamical properties. It also a good starting point for many other more advanced techniques.
During an MD one expects the density to change more or less continously, thus it is possible to use the solutions of the previous steps to create the new initial guess for the density (and wavefunctions for OT).
Indeed cp2k can use different extrapolations techniques (FORCE_EVAL%DFT%QS).
For MD a good guess is PS with an EXTRAPOLATION_ORDER of 3. If you are doing something where from one step to the other there is little continuity, (geometry optimization, montecarlo), then something with little continuity like LINEAR_PS or just the previous density (USE_PREV_P) might be better. If you change the particle number then a USE_GUESS that restarts from scratch with the restart method choosen in FORCE_EVAL%DFT%SCF%SCF_GUESS. The ASPC extrapolation is used in connection with langevin and a reduced SCF convergence (just on SCF step) and hystory restart to simulate big systems.
With the printkey MOTION%PRINT%TRAJECTORY you can control the output of the trajectory. The name of the file is by default projectname-pos-1.xyz and projectname is whatever you have written in GLOBAL%PROJECT. In the prinkey you can also chenge the format from xyz to something else.
An interesting file to check is the projectname-1.ener file, in it you can find the following columns:
md_step, time[fs], e_kin [hartree], temp[K], e_pot [hartree], e_tot [hartree], cpu_time_per_step [s]
You should always check it and look at how the system equilibrates.
The .ener file and other interesting trajectory files are all controlled with subsections of MOTION%PRINT.
If the MD has to be trusted then one should be sure that the trajectory can be trusted. Actually a simple, and very sensitive test that there are no big technical errors is to perform an NVE trajectory and look at the energy conservation. The enrgy conservation has normally two features, short time oscillations (that are larger when the system is not yet equilibrated) and a long range drift. If you express these in K/atom, then you can compare them with the temperature that you are simulating at. Another (equivalent) possibility is to express them as fraction of the kinetic energy. The oscillation and drift (normally per ps, but it also depends on how many picoseconds you want to simulate, and if you want an NVE trajectory) should be small with respect to the kinetic energy (1% or less is a good value, but obviously it depends on the accuracy that you want to achieve, more might be acceptable, or less needed).
To improve the energy conservation one can either improve the forces with EPS_SCF and EPS_DEFAULT, improve ñ increasing the cutoff, and reduce the timestep. For the timestep a good value is normally around 1 fs, but if you are simulating high temperature or special system you might need to reduce it. Normally having good forces and larger timestep is advantageous. Obviously the timestep cannot be smaller than the frequency of the highest oscillations of the system (typically H atoms, that is the reason sometime D atoms are used).
Conserving the total energy well is not enough if the system is very heterogeneous and the interesting part is small. In that case there is the risk that even a large error on it might pass unnoticed. If this is to be excepted (for example the atom with the sharpest gaussian basis set is in the interesting part) then checking that part of the system (oscillations, kinetic energy,...) is advisable.
For the result to be meaningful the system should be equilibrated, and not in a very unlikely state. If one is doing shooting or transition path sampling then this is less true, but still the system should not be in a very high-energy state.
So at the beginning you should build your system in some meaningful way, or from classical simulations. To have a better startingpoint you can optimize the energy (if your system is small), you anneal it, but it is not always needed. Then you can equilibrate it at a given temperature.
To equilibrate a system one can use a crude velocity rescale when he is too far away from the goal temperature (as established by MOTION%MD%TEMP_TOL).
Doing an NVE simulation with TEMP_TOL is better way to equilibrate than using the NVT ensamble that uses the Nose-Hoover chain thermostats and might give spurios effects if you are far from equilibrium, as it tryes to conserve some extended quantity.the thermostat can store energy that he will give back at some later point.
It should be taken into account that the temperature is not constant, but does oscillate in the NVE ensamble, these oscillations are connected with the specific heat and are inversely proportional with sqrt(N) where N is the system size.
A faster way to equilibrate the system is to use LANGEVIN as MOTION%MD%ENSEMBLE, and use a GAMMA of 0.001 [1/fs] (or larger if you want to equilibrate faster). Langevin introduces a viscous dissipative force and a random forces that are in equilibrioum at the given temperature. For equilibration purposes (and not to simulate removed degrees of freedom, like a solvent), the smaller gamma the longer it takes to equilibrate, and the closer the trajectory is to an NVE trajectory, the larger gamma the more the "environment thermostat" influences the system. Also With langevin if your system is in a bad initial configuration it is a good idea to set a TEMP_TOL, or first anneal a little the system.
The Cp2k program is optimized for periodic calculations, and those are the default, still it is possible to perform also other kinds of calculations. First an important thing to note is that by convention the cell goes from 0 to h, not from -h/2 to h/2, thus if ypou want to put a cluster at the center of the cell you need to put it in h/2 not in 0.
The periodicity of the neighboring list, and thus of the gaussian collocation/integration is controlled by FORCE_EVAL%SUBSYS%CELL%PERIODIC.
The poisson solver (electrostatics) is controlled in the section
FORCE_EVAL%DFT%POISSON. The possible poisson solver are choosen with the POISSON_SOLVER keyword, that can have the values
PERIODIC, which can do only fully periodic solvers, ANALYTIC which can do 0, 1d and 2d periodic solutions using analytical green functions in the g space (slow convergence), MT Martyna Tuckermann decoupling, that interacts only with the nearest neighbohr, beware these are completely wrong in the cell is smaller than twice the cluster size (with electronic density),MULTIPOLE, uses a schema that fits the total charge with one gaussian per atom[Blo95].
FORCE_EVAL%DFT%POISSON%PERIODIC, and the poisson solver choosen should support it.
@Article{ lippert:1997,
author = {G. Lippert and J. Hutter and M. Parrinello},
title = {A hybrid {Gaussian} and plane wave density functional
scheme},
journal = mp,
volume = {92},
number = {3},
pages = {477--487},
year = {1997}
}
[VKM04]
@Article{ qs_review:2004,
author = {Joost VandeVondele and Matthias Krack and Fawzi Mohamed
and Michele Parrinello and Thomas Chassaing and J{\"u}rg
Hutter},
journal = cpc,
title = {{\sc Quickstep}: fast and accurate density functional
calculations using a mixed Gaussian and plane waves
approach},
pages = {103--128},
volume = 167,
year = {2005}
}
[VH03]
@Article{ vandevondele:2003,
author = {J. VandeVondele and J. Hutter},
title = {An efficient orbital transformation method for electronic
structure calculations},
journal = jcp,
volume = {118},
number = {10},
pages = {4365--4369},
year = {2003}
}
[Blo95]
@Article{bloechl:1995,
author={P. E. Bl{\"o}chl},
journal=jcp,
year=1995,
volume=103,
issue=17,
title = {Electrostatic decoupling of periodic images of plane-wave-expanded
densities and derived atomic point charges},
pages={7422--7428}
}
[GTH96]
@Article{ goedecker:1996,
author = {S. Goedecker and M. Teter and J. Hutter},
title = {Separable dual-space {Gaussian} pseudopotentials},
journal = prb,
volume = {54},
number = {3},
pages = {1703--1710},
year = {1996}
}