This is an old revision of the document!
Table of Contents
Ab initio molecular dynamics
Once a standard GGA simulation has been set up, doing ab initio MD is easy. Here we prepare a simulation of bulk liquid water, a system that has been studied a lot with CP2K (e.g. 10.1063/1.1828433 or 10.1021/jp901990u). The second illustrates convincingly why dispersion corrections are essential.
The first goal is to setup a simulation in production mode by reducing output, and enabling restarting. The second goal to understand the produced .ener file and do some basic analysis of the trajectory with VMD.
AIMD of water
Topics:
- MD section (timestep)
- Thermostat (NVE, NVT, NPT)
1st task: prepapre inputs for MD
Start from the mode1.inp
input file of the previous exercise, rename the file to water.inp
. You'll also need a new file HFX_BASIS
from the cp2k/data
directory.
Change keywords to :
PROJECT WATER
RUN_TYPE MD
BASIS_SET_FILE_NAME HFX_BASIS
MINIMIZER DIIS
ABC [angstrom] 12.42 12.42 12.42
COORD_FILE_NAME water.xyz
BASIS_SET DZVP-GTH
Insert the following blocks in their respective sections:
Reduce output, and limit job time to 5min.
&GLOBAL ! limit the runs to 5min WALLTIME 300 ! reduce the amount of IO IOLEVEL LOW &END GLOBAL
! do not write the wfn restart file every step, for large systems this is slow
&SCF ! do not store the wfn during MD &PRINT &RESTART OFF &END &END &END SCF
! tune the output, switch off unneeded keywords ! only sample output of others.
&PRINT ! at the end of the SCF procedure generate cube files of the density &E_DENSITY_CUBE OFF &END E_DENSITY_CUBE ! compute eigenvalues and homo-lumo gap each 10nd MD step &MO_CUBES ! compute 4 unoccupied orbital energies NLUMO 4 NHOMO 4 ! but don't write the cube files WRITE_CUBE .FALSE. ! do this every 10th MD step. &EACH MD 10 &END &END &END
make sure trajectories and restart files are dumped as needed.
&MOTION &PRINT &TRAJECTORY &EACH MD 1 &END EACH &END TRAJECTORY &VELOCITIES OFF &END VELOCITIES &FORCES OFF &END FORCES &RESTART_HISTORY &EACH MD 500 &END EACH &END RESTART_HISTORY &RESTART BACKUP_COPIES 3 &EACH MD 1 &END EACH &END RESTART &END PRINT &END MOTION
Now, run the input and check for successful completion of the job (look for the timing report or the ENDED AT
line in the output).
During MD CP2K has generated various files named WATER-1.restart
, WATER-pos-1.xyz
, WATER-1.ener
. These are a restart file for the MD (open this file in an editor to see that this is actually a regular input file), the trajectory of the MD, and a file summarizing the energies (potential, kinetic and conserved quantity). Before we analyze these files, we'll resubmit a job for the continuation of the MD.
modify the input to enable restarting:
&EXT_RESTART RESTART_FILE_NAME WATER-1.restart &END
If time permits, increase the WALLTIME (and the corresponding time limit in the job script), and run a longer job.
bsub
, or –dependency of sbatch
Analysis
While running the MD simulations, it is useful to check that all is fine. Here, we do some basic analysis, to look at the structure and dynamics of the liquid.
2nd task: visualize the .ener file
A first quick check can be performed using the file WATER-1.ener
which contains basic information
# Step Nr. Time[fs] Kin.[a.u.] Temp[K] Pot.[a.u.] Cons Qty[a.u.] UsedTime[s] 0 0.000000 0.273612846 300.000000000 -1102.629448594 -1102.355835748 0.000000000 1 0.500000 0.279633819 306.601634956 -1102.634728486 -1102.356032711 70.082937483 2 1.000000 0.278176228 305.003473822 -1102.643688340 -1102.356285321 11.608515253 3 1.500000 0.280393422 307.434493169 -1102.653080703 -1102.356547289 11.607597935 4 2.000000 0.282889483 310.171274373 -1102.655862600 -1102.356593452 11.617385623 5 2.500000 0.294372846 322.762089451 -1102.653391721 -1102.356505399 11.665471402
This can be easily visualized with gnuplot, for example for the conserved quantity (plotting the second vs the sixth column) :
gnuplot> plot './WATER-1.ener' u 2:6 w lp
To judge if this is actually well conserved, compare to the potential energy:
gnuplot> plot './WATER-1.ener' u 2:5 w lp gnuplot> replot './WATER-1.ener' u 2:6 w lp
If the constant of motion is not well conserved, try to
- Make
EPS_SCF
tighter (to reduce drift) - Make
TIMESTEP
shorter (to reduce fluctuations) - Play with
EXTRAPOLATION_ORDER
(to reduce drift and or instabilities)
To judge if a system is well equilibrated is not easy. At least the temperature and the potential energy of the system must oscillate around an average and be free of long term drift. As a rule of thumb, discard 1/3 of the trajectory, use 2/3 for data analysis.
3rd task: visualize/analyze the trajectory file
We will use vmd to analyze the trajectory file. Note that the generated trajectory is only a few 100s of fs, typically, 10s of ps are needed for even for basic properties.
Start vmd
vmd WATER-pos-1.xyz
g(r)
In the menu go to :
Extensions/Analysis/Radial Pair Distribution Function g® Utilities/Set unit cell size dimensions
1st, set the unit cell as needed. Now improve the Graphics/Representations to also show neighboring unit cells and visualize hydrogen bonds.
2nd, compute the O-O pair distribution function (Selections=name O
) and similar for the O-H pair distribution function (including their integrals).
How many neighbors does a given water molecule have on average (3, 3-4, 4, 4-5, 5)?
IR spectrum
Based on the time evolution of the dipole of the system, the IR spectral density can be estimated. To estimate the dipole from AIMD, wannier centers need to be computed. This is out of scope of the current tutorial (TODO: find link). We employ a simple approximation, namely classical point charges for the water molecules. In this context the approximation is reasonable.
Create the following file
- charges.dat
O -1.2 H +0.6
Go to Extensions/Analysis/Spectral density calculator. Select the proper molecule (WATER-pos-1.xyz), adjust the timestep (0.5fs), and the maximum frequency (6000 cm^-1). Utilities/Load name↔charge map from file. Compute spectrum.
Where do you expect the OH stretch to be ? Is this reproduced ?
4th task: simple ions in solution
Introduce an ion in your system, and equilibrate this system. Study its dynamics and solvation structure.
The easiest way to do so is to replace one or more water molecules (depending on the size of the ion) by the ion in question. Obviously, the configuration produced in this way is far from equilibrium, and must be run for a while before it is representative.
Entertaining is to turn one H2O in H+, do you see Eigen and Zundel states and the Grotthuss mechanism ?
Required files
- water.xyz
192 water with unit cell: ABC [angstrom] 12.42 12.42 12.42 O 3.8585873763 -4.7533175213 5.5091974759 H 4.2109950951 -5.5568705259 5.9994585948 H 3.0947084560 -4.3172663108 5.9135950646 O -3.5460761891 1.7456237131 4.2880212708 H -3.3572582677 0.8369944362 4.1529456220 H -3.6089143450 1.9436976972 5.2698392012 O 4.5691359319 2.9343274990 -7.0999718127 H 3.6785962350 3.2291072409 -6.8049831901 H 4.4953665356 2.7845914189 -8.0458027314 O -2.5731801314 -0.1063102257 -3.4147140374 H -2.0386687011 0.4155234251 -2.8289645480 H -1.8891606976 -0.6148201666 -3.9280631618 O -3.7625034284 -7.7105940327 5.6279977389 H -4.6343536576 -7.4199803039 5.9874790347 H -3.8520633915 -7.6348492192 4.6607806353 O -6.7131259426 -0.0759204976 2.3176844292 H -6.0415850051 -0.4560400697 3.0079983076 H -6.2582643134 0.7343068494 1.9582947358 O 4.0109535867 -5.5817692668 -0.4126912859 H 3.3321311866 -5.6166906662 0.3076935926 H 3.9051538355 -4.7533442556 -0.8369521875 O -2.6692626967 -3.0810510613 -1.2259349448 H -3.2486322997 -2.4491708867 -1.7386930105 H -2.1918958173 -2.5988518360 -0.5763067756 O 5.4576244064 -3.1688044689 3.1002596795 H 4.6854987571 -3.2457995297 3.6224179934 H 5.2579372917 -3.0507495002 2.1311396222 O 2.0078792352 3.8025888294 -6.6441198936 H 1.3809036112 3.0702729257 -6.8038689306 H 1.8366067800 4.0666291804 -5.7099975874 O 4.0608830110 -3.8929207302 -3.0596772434 H 3.2332712864 -3.7206669725 -3.6953758544 H 4.3366899606 -2.9825098616 -2.9848217710 O 1.9773593054 -1.8442291306 1.0225667179 H 1.5650644885 -2.4713604568 1.6495937537 H 1.6194813185 -1.9133123969 0.0656946095 O 0.7338447631 -5.7432877136 -5.0100685957 H 1.4795768454 -6.3299226862 -4.8189325814 H 0.4491149358 -5.7095685669 -5.9473559749 O -7.5166459589 -8.6781062088 -3.2371569882 H -7.8737663294 -9.6027219294 -3.0814153209 H -7.1323941392 -8.5788960418 -2.3467663280 O 3.0089956358 -1.6206934676 5.1704516391 H 3.8547864750 -1.8022862145 5.6866350548 H 3.1700832113 -0.8185381416 4.6102432456 O -0.0943492744 1.8634710749 -6.6861347029 H 0.2389831969 1.0552334471 -6.3170377387 H -0.0534604238 1.7370469556 -7.6527831571 O 1.3803001006 -3.0815472275 3.6598466091 H 2.0052752091 -2.6793618758 4.3169228939 H 1.0973942250 -3.9666380514 3.9813401483 O -2.3279129364 2.8060791527 2.2347417439 H -1.6230445751 2.1315230876 2.1866163236 H -2.7609234795 2.5819568350 3.0889511952 O 0.8745131234 -0.6611187482 -5.6004046460 H 1.7204015713 -0.8189656405 -6.1098612966 H 1.1964117837 -0.2405057963 -4.7549061774 O -1.1866462102 -0.7164030048 0.2437849238 H -1.1340006348 0.2325710411 -0.0726156939 H -0.9335438317 -0.7854778370 1.2011335414 O -3.8861935106 -9.0129285534 -3.1318571683 H -4.4983492231 -8.2801255892 -3.2160220896 H -2.9183857578 -8.7301395642 -3.1224505015 O 4.9784414381 4.5574106629 2.1605978288 H 4.5742506248 4.0690251974 1.4375040342 H 4.2106782789 4.9079438424 2.6736093641 O 4.2074812608 1.0294752439 -3.2134938111 H 4.6250866171 0.8524308220 -4.0648666086 H 4.4221519308 0.1920392568 -2.7364360799 O -5.2106345443 -6.3463699492 -2.6806456410 H -5.3801630197 -6.0861826108 -1.7737611737 H -5.9985546753 -6.0894306119 -3.1521879087 O 4.8272006144 -2.1493089974 0.6957555837 H 5.1698114049 -1.3843105780 1.1680731517 H 3.8731403616 -2.0152767607 0.7901517878 O 2.1084974235 0.9406557281 1.6964699575 H 2.7635704789 1.5106519967 1.2582480434 H 2.3138417141 0.0791504768 1.2943222699 O -5.7579045611 -5.0935991816 -0.3613728277 H -6.6922016300 -5.4438912590 -0.2434456391 H -5.9339541154 -4.1852263482 -0.6893230982 O -4.8501183696 4.6347512586 2.6323657871 H -5.8186923573 4.7291496971 2.3648078237 H -4.7686962808 3.6685062675 2.2489112658 O 5.5123889390 0.8523825517 -5.5927914459 H 6.4009274201 1.1961666201 -5.3454982775 H 5.0326268549 1.4571097750 -6.2820870000 O 0.7365184732 -1.9176950430 -1.2724577473 H 0.4366627686 -2.6685959517 -1.8452499204 H -0.0280723827 -1.6346960324 -0.7323636982 O 3.9005939808 2.6221972296 0.0824981180 H 4.6603628952 2.4136405159 -0.4885023274 H 3.3619307241 3.2465911837 -0.4599512808 O -1.0991624885 -1.3693412589 2.9498880034 H -0.3354058556 -2.0369012425 3.1075024163 H -1.8473458549 -1.9916364646 3.2032143743 O 2.0718451094 -5.6501678174 1.3629079156 H 1.2043618452 -6.1369822674 1.2560686630 H 2.4507495844 -5.8859092245 2.2558212721 O -0.6280310316 -3.8273493362 -2.8460859716 H -1.5583639221 -3.7769636412 -2.5811140145 H -0.5112292363 -4.5879523025 -3.3982225959 O -5.2918433010 2.0621670074 1.4353011132 H -4.4549447831 1.8776050679 1.9014549175 H -5.1335936869 1.9198214761 0.4197228743 O 1.9599719653 -3.4245632717 -5.2673604360 H 1.3225562775 -4.2123092010 -5.0388832852 H 1.2999292168 -2.7711468268 -5.4986783863 O -1.1933205846 3.9368513830 -3.3737465077 H -1.2003012694 3.9803649971 -4.3582935942 H -0.7402505311 4.6500956055 -2.9405594580 O -5.4243979802 -5.0029807582 4.2784625506 H -6.1164034705 -4.5127003053 3.7817315683 H -4.7403199519 -5.3498607764 3.6920944114 O 2.9732984993 5.7937546317 3.6041830736 H 3.4176710045 6.4726708296 4.1695014822 H 2.4485201643 5.2041088060 4.2250869218 O -4.0511295134 1.5445925671 -5.0626498585 H -3.5387103609 0.7974835831 -4.6210546124 H -4.0333821170 2.3518494025 -4.4301366268 O 3.4847726839 0.5997692637 4.0091486817 H 2.7896502362 0.5375110591 3.2686671028 H 4.3189983631 0.4153324309 3.4395276451 O -4.6207476869 -1.0117618538 4.0048952776 H -5.3477593066 -1.2443059233 4.6253237945 H -4.2163881176 -1.8844792105 3.8428058821 O -0.2165387154 1.2751467697 3.1167167744 H 0.6371782887 1.4136801810 2.7286969429 H -0.3894124526 0.3225841245 3.0508790291 O -1.9915351632 5.0697083632 0.6000194332 H -2.3376502958 4.2340876605 1.0013898546 H -2.6520248188 5.7471163392 0.8506616742 O -1.3075474036 -1.8956122645 -4.9542920356 H -0.8701251561 -2.4028189929 -4.2563741889 H -0.5949693850 -1.3532124979 -5.4532434800 O -5.1373452268 1.8062068789 -1.2867681226 H -4.6981592560 2.4597481938 -1.8505668997 H -5.4245526684 1.1256927260 -1.8633945543 O -7.3836561358 -1.2119873462 -1.9506599575 H -6.4212881329 -1.0371432188 -2.1463799566 H -7.4808782211 -1.5310477918 -1.0339482458 O -0.4639445925 1.4048471789 -1.7228637494 H -0.7867348868 2.1866499220 -2.1962304023 H 0.2743927320 0.9864157199 -2.2215845685 O -0.0298107043 -6.5210921839 -1.0933643185 H 0.0021424143 -5.5311453888 -1.1226646288 H -0.8031132631 -6.7435142146 -0.4648562320 O 5.7090935537 -1.6927086580 6.0185486672 H 6.0227718591 -2.5102376528 6.4195407580 H 5.6793497711 -0.9012395738 6.6521961143 O -2.5962982991 -5.7683972037 -4.2195771395 H -2.6049028140 -6.7077669400 -4.4728336739 H -3.4354537463 -5.5645841736 -3.8483514813 O -4.9449687319 -1.4086720170 -3.2473392351 H -4.8913021653 -2.2217322989 -3.7674224672 H -3.9689548003 -1.0964857745 -3.2341980948 O 2.2983754229 -8.0809494482 -1.2132412662 H 1.3003681505 -8.0636283692 -1.1143923954 H 2.5161940132 -7.1053673603 -1.1768519202 O 6.7754357543 -3.8965331032 -5.3359261145 H 7.5499418595 -4.1003001805 -5.8213622687 H 6.0430158524 -4.4328146397 -5.6479610986 O -0.3424844631 4.9658112476 2.9508285728 H -0.6418240309 4.1672557959 3.4716379980 H -0.6759740908 4.8699377537 2.0551376194 O -2.8056288807 -4.0048234272 6.2239057927 H -2.3850599469 -4.6322887909 6.8455564859 H -2.3797353004 -3.1498233269 6.4330163166 O -2.7925113864 -3.4732237158 3.4612331596 H -3.1129543056 -3.8620641144 4.2746132111 H -3.2856093678 -4.0271447136 2.8274407488 O 2.4867711476 4.6803618907 -3.9511208166 H 3.3971373036 4.5081635205 -3.6298163906 H 1.9790712168 4.6240682409 -3.1152437525 O 1.6590765104 0.0523137912 -3.1069189687 H 1.5113564742 -0.5996769003 -2.3563690233 H 2.5068721250 0.5264194699 -3.0661523880 O -1.2023048654 -8.1603479212 -6.3670700282 H -0.9055157967 -9.0690551157 -6.5705273597 H -2.1798096138 -8.0740016376 -6.5464905196 O -0.4761071273 8.3851620636 0.5987466087 H 0.3632441731 8.3096547993 1.1058350782 H -1.2852386123 8.0753854402 1.0587255936 O 5.9065449362 -6.7989771803 -6.6004398209 H 5.4130154976 -7.6233171741 -6.9137876216 H 6.2177506878 -6.2344541440 -7.3198720135 O -0.3594330297 -5.2622440761 4.6704297767 H -0.4156127476 -6.0988819305 4.1367481919 H -1.2931127588 -5.1547968831 4.9497584977 O -3.8273131578 -5.5279235466 1.7633664403 H -4.1388543582 -6.3621998573 2.1513804320 H -4.3955669691 -5.3216936959 1.0405948547