<h1>CP2K Tutorial</h1>
<p>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 <a href="mailto:fawzi@gmx.ch">fawzi@gmx.ch</a>. If you have further questions about cp2k, then go to <a href="http://cp2k.berlios.de/">http://cp2k.berlios.de/</a> 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.</p>
<h2>Starting out/Installing</h2>

<h3>Berlios</h3>
One of the main starting points is the Berlios site of cp2k <a href="http://cp2k.berlios.de/">http://cp2k.berlios.de/</a>. On it there is information on cp2k, documentation, installation, regtests, mailinglists, and you can contact the developers.

<h3>Getting and installing CP2K</h3>

<p>One of the easiest ways to get cp2k if one has cvs installed (and it is a good idea) is <br>
<code>touch ~/.cvspass<br>
cvs -d :pserver:anonymous@cvs.cp2k.berlios.de:/cvsroot/cp2k login<br>
cvs -d :pserver:anonymous@cvs.cp2k.berlios.de:/cvsroot/cp2k co cp2k<br>
</code>
You can add <code> -D YYYY-MM-DDTHH:MM</code> to specify the date.</p>

<p>To compile it you have to adapt the architecture file in <code>cp2k/arch</code>. Executing <code>cp2k/tools/get_arch_name</code> you can see your architecture, changing the environment variable <code>FORT_C_NAME</code> you can change the compiler on some architectures.</p>
<p>As <code>make</code> you should use the gnu make, and the architecture can be changed also with <code>make ARCH=my_arch_name</code>.</p>
<p>The kind of build (sdbg, sopt, pdbg, popt) can be changed with <code>VERSION=build_type</code> (actually using the build type as target also works).
Targets can be: build, lib, clean, realclean, distclean.</p>

<h3>Regtests</h3>
<p>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.</p>
<p>To setup the regtests you have to
    <ul>
    <li>create a test directory with a copy of the <code>cp2k</code> directory and of <code>cp2k/tools/do_regtest</code> into it for example with<br>
        <code>mkdir tst<br>
            cd tst<br>
            cvs -d :pserver:anonymous@cvs.cp2k.berlios.de:/cvsroot/cp2k co cp2k<br>
            cp cp2k/tools/do_regtest .
        </code>
    </li>
    <li>
        edit the do_regtest script, you will want to set up the things so that the compiler (and maybe the mpi) are available, fix FORT_C_NAME, the base directory (the one that contains <code>cp2k</code> and <code>do_regtest</code>), decide the version you want to check and set dir_triplet (normally to what is returned by <code>cp2k/tools/get_arch_name</code>).
        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 <code>datum</code> like
<code>date -u '+%Y-%m-%dT%H:%M:%S'</code> to avoid empty spaces in it.
    </li>
</ul>

<h3>Input Documentation</h3>

<h4>Documentation</h4>
<p>Documentation of cp2k can be found on the website <a href="http://cp2k.berlios.de/">http://cp2k.berlios.de/</a> and can also be generated by the cp2k executable itself using the flag <code>--html</code> or <code>--xml</code> (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.</p>

<h4>Input Syntax</h4>
<ul>
<li>the input is subdivided in sections that start with "&amp;" and end with "&amp;END"</li>
<li>each section can contain other sections and keywords followed by values</li>
<li>both the section names and the keyword names are case insensitive</li>
<li>the special characters "#" and "!" denote the beginning of a comment that goes until the end of the line</li>
<li>If the value of the keyword is experssed in units, perceding the value by the units in square bracket changes them. In the square brackets the units should have no space. For example "[m]" switches to meters.
<li>A long line of input can be split using the "\" character at the end of the line</li>
<li>keywords with a logical value (flags), often have a value that is false when the keyword is missing; the value true when one writes just the keyword (lone_keyword_value), or they can be given an explicit value after the keyword as normal (for example <code>KEYWORD FALSE</code>).</li>
<li>Sections can have a default parameter that comes right after the section value. For example the <code>OT</code> section has a logical default parameter that controls the activation of the OT method. So,
<br><code>&amp;OT
<br>...
<br>&amp;END OT</code><br>
activates the OT method and
<br><code>&amp;OT FALSE
<br>...
<br>&amp;END OT</code><br>
(or the absence of the OT section) deactivates it.</li>
</ul>

<h4>Input Structure</h4>
<p>The 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 <code>FORCE_EVAL</code> section.
If you use the ab-initio DFT part then the <code>DFT</code> and <code>QS</code> subsections will probabily be of interest for you.
The method to use is choosen with <code>FORCE_EVAL%METHOD</code> keyword, i.e. the keyword <code>METHOD</code> in the <code>FORCE_EVAL</code> section. So you can easily switch from one method to the another if you gave meaningful parameters to both.
</p>

<p>Then parameters for the kind of calculation (run type) to do using that method (energy,force,molecular dynamics,langevin,...) are defined in the <code>MOTION</code> section</p>

<p>Finally in the <code>GLOBAL</code> section you can choose the run type to perform with the keyword <code>RUN_TYPE</code>, the name of the project (<code>PROJECT</code>), used to generate most of the output files, the <code>PRINT_LEVEL</code> keyword that controls the general level of output, and other general flags. The <code>GLOBAL%PRINT%PHYSCON</code> keyword enables a printout of the physical constants used by cp2k.</p>

<h4>Print Keys</h4>
<p>To control the output of information or properties often <i>printkeys</i> 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 <code>PRINT</code> section.<p>

<p>To understand how to control the output it is necessary to first understand the <i>iteration level</i>. 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.</p>
<p>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</p>

<table border="1">
<tr><td>section parameter</td>
<td>	Level starting at which this property is printed, valid values are
<code>ON, OFF, SILENT, LOW, MEDIUM, HIGH, DEBUG</code>
	No section means the default value, the presence of the prind sections implies a value of <code>ON</code> (always printed), if one wants to selectively switch off the printkey it is possible the give it the explicit value OFF, i.e. <code>&amp;SECTION OFF</code>
</td>
</tr>
<tr>
<td>EACH</td>
<td>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)</td>
</tr>
<tr>
<td>ADD_LAST</td>
	<td>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</td>
</tr>
<tr>
<td>COMMON_ITERATION_LEVELS</td>
<td>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)
</td>
</tr>
<tr>
<td>FILENAME</td>
	<td>controls part of the filename for output.
<ul><li>Use "__STD_OUT__" (exactly as written here) for the screen or standard logger.</li>
<li>Use "filename" to get "projectname-filename[-middlename]-iteration.extension"</li>
<li>Use "./filename" (or any path that contains at least one "/") to get "./filename[-middlename]-iteration.extension"</li>
<li>Use "=filename" to get exactly filename as output file (please note that this can be dangerous because different output/iterations might write on the same file</li>
</td>
</tr>
</table>

<h4>Restarting</h4>
<p>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.</p>

<p>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 <code>EXT_RESTART</code> giving the name of the external input with the <code>EXTERNAL_FILE_NAME</code> 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 <code>RESTART_VEL FALSE</code>.</p>

<p>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 <code>FORCE_EVAL%DFT%RESTART_FILE_NAME</code>.</p>

<h3>Sample Inputs</h3>

<p>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.</p>

<h2>GPW</h2>
<p>The Gaussian Plane Waves method [<a href="#Lip97">Lip97</a>] 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.</p>
<p>In GPW the whole density is transferred to plane waves, and one has the density <i>n(r)=sum<sub>i j</sub>P<sup>i j</sup> phi<sub>i</sub>(r)phi<sub>j</sub>(r)</i> in the gaussian basis set and the density <i>&ntilde;</i> taken on a uniform grid.</p>

<p>Then <i>&ntilde;</i> is used to caluclate the hartree potential <i>V<sub>H</sub></i> (via an FFT based poisson solver) and the exchange and correlation potential <i>V<sub>xc</sub></i>. These potential are then transferred back to gaussian basis set by integrating them with <i>phi<sub>i</sub>(r)phi<sub>j</sub>(r)</i>. To make the collocation and integration perfectly consistent with each other one can set <code>FORCE_EVAL%DFT%QS%MAP_CONSISTENT</code> this normally adds only a very small overhead to the calculation and is </p>

<h3>Cutoff</h3>
<p><i>n</i> and <i>&ntilde;</i> are not equal, and this introduces an error in the calculation. <i>&ntilde;</i> converges toward <i>n</i> 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).</p>

<p>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 <i>0.5(pi/dr)**2</i> where dr is the gridspacing in Bohr. The characteristic length of a gaussian with exponent A is given by <i>1/sqrt(a)</i> (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 <i>alpha<sub>H</sub>=47.7</i> a good cutoff for doing MD is 280 Ry.</p>

<p>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 [<a href="#GTH96">GTH96</a>], these are harder than other pseudopotentials, but also more accurate.</p>

<h3>Smoothing</h3>
<p><i>&ntilde;</i> 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 <code>XC_SMOOTH_RHO NN10</code> and <code>XC_DERIV SPLINE2_SMOOTH</code> (in the <code>FORCE_EVAL%DFT%SC%XC_GRID</code> 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.</p>
<p>Methods that do not redefine the total energy are <code>XC_SMOOTH_RHO NONE</code> and <code>XC_DERIV</code> equal to either <code>PW, SPLINE3</code> or <code>SPLINE2</code>. Thes are listed from the one that assumes more regularity (<code>PW</code> the the one that assumes less regularity <code>SPLINE2</code>. Normally <code>SPLINE2</code> is a good choice, but for high cutoffs (600 Ry for water) <code>SPLINE3</code> is better. The default (<code>PW</code>) is not bad, but generally inferior to the others.</p>

<h3>Basisset</h3>
<p>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.</p>
<p>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.</p>
<p>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.</p>
<p>The basis set to use for an element is set in the <code>FORCE_EVAL%SUBSYS%KIND</code> section with its name with the <code>BASIS_SET</code> keyword.</p>
<h3>Pseudo Potential</h3>
<p>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.</p>
<p>The pseudopotential to use for an element is set in the <code>FORCE_EVAL%SUBSYS%KIND</code> section with its name.</p>
<h3>XC selection</h3>
<p>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.</p>
<p>In cp2k the xc functional is selected in the <code>FORCE_EVAL%DFT%TDDFPT%XC</code> section. Common functional combinations can be selected directly as section parameters of the <code>XC_FUNCTIONAL</code> for example with the <code>POTENTIAL</code> keyword.</p>
<br><code>&amp;XC_FUNCTIONAL BLIP
<br>&amp;END XC_FUNCTIONAL</code><br>
but more complex combination have to be choosen by directly setting the underlying subsections.</p>

<h3>WF Convergence</h3>
<p>Apart from cutoff and basis set with the GPW method there are two other importan switches that control the convergence of the method. One is <code>FORCE_EVAL%DFT%QS%EPS_DEFAULT</code> that controls the convergence of integral calculation, overlap cutoff, neighboring lists... and sets a slew of other <code>EPS_*</code> keywords so that the KS matrix is calculated with the requested precision. The other is <code>EPS_SCF</code> which controls the convergence of the scf procedure. In general one should have an error <i>e<sub>f</sub></i> on the forces one should have <i>sqrt(eps_scf) &lt; e<sub>f</sub></i> and <i>sqrt(eps_default) &lt; eps_scf</i>, and in general around 1.e-3 is the error that you need to have reasonable MD.</p>

<h3>ot / diag</h3>
<p>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.</p>
<p><code>FORCE_EVAL%DFT%SCF%SCF_GUESS RESTART</code> is a parameter that controls how the initial density is built, normally one uses either <code>ATOMIC</code> or <code>RESTART</code>. During an MD this setting normally influences only the initial startup.</p>
<p>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 <code>FORCE_EVAL%DFT%SCF%MIXING</code> factor. If the new density is close enough to the old density <code>FORCE_EVAL%DFT%SCF%EPS_DIIS</code> then the DIIS procedure is activated.</p>
<p>Diagonalization works well, but for difficult or big systems the OT method [<a href="VH03">VH03</a>] 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 <code>FORCE_EVAL%DFT%SCF%OT</code> is enough.</p>
<p>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.</p>
<p>The speed of convergence of the OT method the preconditioner choice (<code>FORCE_EVAL%DFT%SCF%PRECONDITIONER</p>) is important. <code>FULL_KINETIC</code> is a good choice for very large systems, but for smaller or difficult systems <code>FULL_ALL</code> is normally better. The quality of the preconditioner is dependent on how close one is to the solution. Especially with expensive preconditioners like <code>FULL_ALL</code> 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 <code>FORCE_EVAL%DFT%SCF%MAX_SCF</code> to a small value, and activating the outer scf loop addint the <code>FORCE_EVAL%DFT%SCF%OUTER_SCF</code> section, and setting its <code>EPS_SCF</code> equal to the one of the scf loop, and chooing a outer <code>MAX_SCF</code> big enough (the total number of iteration is the product of internal and external MAX_SCF.</p>
<p>Also the minimizer of the OT method can be changed with<code>FORCE_EVAL%DFT%SCF%OT%MINIMIZER</code>, the default is <code>CG</code> that uses a very robust conjugated gradient method. Less roubust but often faste is the <code>DIIS</code> method.</p>

<h2>MD</h2>

<p>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.</p>
<p>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 (<code>FORCE_EVAL%DFT%QS</code>).</p>
<p>For MD a good guess is <code>PS</code> with an <code>EXTRAPOLATION_ORDER</code> 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 <code>LINEAR_PS</code> or just the previous density (<code>USE_PREV_P</code>) might be better. If you change the particle number then a <code>USE_GUESS</code> that restarts from scratch with the restart method choosen in <code>FORCE_EVAL%DFT%SCF%SCF_GUESS</code>. The <code>ASPC</code> extrapolation is used in connection with langevin and a reduced SCF convergence (just on SCF step) and hystory restart to simulate big systems.</p>

<h3>Trajectory</h3>
<p>With the printkey <code>MOTION%PRINT%TRAJECTORY</code> you can control the output of the trajectory. The name of the file is by default <b>projectname-pos-1.xyz</b> and projectname is whatever you have written in <code>GLOBAL%PROJECT</code>. In the prinkey you can also chenge the format from xyz to something else.</p>
<p>An interesting file to check is the <b>projectname-1.ener</b> file, in it you can find the following columns:
<br><b>md_step, time[fs], e_kin [hartree], temp[K], e_pot [hartree], e_tot [hartree], cpu_time_per_step [s]</b><br>
You should always check it and look at how the system equilibrates.</p>
<p>The .ener file and other interesting trajectory files are all controlled with subsections of <code>MOTION%PRINT</code>.</p>

<h3>MD Convergence</h3>
<p>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 <i>K/atom</i>, 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 <i>ps</i>, 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).</p>
<p>To improve the energy conservation one can either improve the forces with EPS_SCF and EPS_DEFAULT, improve <i>&ntilde;</i> 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).</p>
<p>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.</p>
<h4>Equilibration</h4>
<p>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. </p>
<p>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.</p>
<p>To equilibrate a system one can use a crude velocity rescale when he is too far away from the goal temperature (as established by <code>MOTION%MD%TEMP_TOL</code>).
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.</p>
<p>A faster way to equilibrate the system is to use <code>LANGEVIN</code> as <code>MOTION%MD%ENSEMBLE</code>, and use a <code>GAMMA</code> 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.
</p>
<h2>gapw</h2>

<h2>Periodicity</h2>
<p>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.</p>
<p>The periodicity of the neighboring list, and thus of the gaussian collocation/integration is controlled by <code>FORCE_EVAL%SUBSYS%CELL%PERIODIC</code>.</p>
<p>
The poisson solver (electrostatics) is controlled in the section
<code>FORCE_EVAL%DFT%POISSON</code>. The possible poisson solver are choosen with the <code>POISSON_SOLVER</code> keyword, that can have the values 
<ul>
<li><code>PERIODIC</code>, which can do only fully periodic solvers, </li>
<li><code>ANALYTIC</code> which can do 0, 1d and 2d periodic solutions using analytical green functions in the g space (slow convergence), </li>
<li><code>MT</code> 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),</li>
<li><code>MULTIPOLE</code>, uses a schema that fits the total charge with one gaussian per atom[<a href="Blo95">Blo95</a>].
</li>
</ul>
The periodicity of the poisson solver (electrostatics) is set in <code>FORCE_EVAL%DFT%POISSON%PERIODIC</code>, and the poisson solver choosen should support it.</p>

<h2>QM/MM</h2>

<h2>References</h2>
<a name="Lip97">[Lip97]</a>
<pre>@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}
}</pre>

<a name="VKM04">[VKM04]</a>
<pre>@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}
}</pre>

<a name="VH03">[VH03]</a>
<pre>@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}
}</pre>

<a name="Blo95">[Blo95]</a>
<pre>@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}
}</pre>

<a name="GTH96">[GTH96]</a>
<pre>@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}
}</pre>