Talk:Mod:Hunt Research Group/cpmd water
Car-Parrinello Molecular Dynamics Simulation of Aqueous Solutions
1. Introduction
1.1 Aqueous solutions
In this lab you are going to run CPMD calculations of systems containing water and chloride. The following paragraphs describe why these systems are interesting to study.
Solvation of the chloride anion in water solution has been the subject of a substantial amount of experimental and theoretical work, in view of the central importance of chloride solutions in biochemical and metabolic processes, geology, atmospheric chemistry, and the chemical industry. Indeed, it has been discovered that airborne aqueous sea-salt aerosols serve as a global source of molecular chlorine and bromine, both in the polluted and in the remote lower marine troposphere. A variety of experimental methods have been used to investigate the hydration structure of chloride ions in water, including neutron and X-ray diffraction techniques, X-ray absorption spectroscopy (XAS) and femtosecond mid-infrared nonlinear spectroscopy. Despite these efforts, a complete microscopic model of the hydration shell structure of aqueous Cl- and of its room temperature dynamics is still being debated.
Difficulties in determining the structure of the solvation environment arise because of the comparable magnitudes of the water-Cl- and water-water interaction energies. Established results from the literature, both experimental and theoretical, reveal a rather inhomogeneous picture of the hydrated ion solvation shell structure. Vibrational spectroscopy experiments indicate that a chloride ion is asymmetrically solvated in clusters containing up to five water molecules. This contrasts with older and less direct experiments based on mass spectroscopy and photoelectron spectroscopy from which a symmetrical solvent cage around the anion was inferred. Experimental measurements such as photoelectron spectroscopy, IR spectroscopy and a combined technique of XAS (extended X-ray absorption fine structure, EXAFS and X-ray absorption near-edge structure, XANES) yield coordination numbers that are significantly scattered, ranging from 3 to 8. The first peak of the Cl-O radial distribution function (RDF) is measured to be in the range of 3.10-3.36 A.
Uncertainty in the coordination number is also evinced from theoretical studies, with variations between 5.1 and 8.4 depending on the simulation techniques adopted. Classical molecular dynamics (MD) simulations have been found to depend strongly on the anion-water and water-water potentials used. On the other hand, the use of different parametrized potentials, i.e. with and without a treatment of molecular polarizability, gave different structural properties. Quantum mechanics / molecular mechanics (QM/MM) simulations of Cl− in water indicate that the hydration structure of the an- ion has considerable flexibility in consequence to the weakness of Cl--water interactions. This leads to a competition between solvation of the ion and hydrogen bonding among water molecules. The QM/MM results clearly indicate the importance of QM treatment in order to correctly describe such a delicate balance between Cl--water and water-water interactions.
Reference: Tongraar, A.; T-Thienprasert, J.; Rujirawat, S.; Limpijumnong, S. Phys. Chem. Chem. Phys. 2010, 12, 10876.
1.2 Car-Parrinello/Born-Oppenheimer molecular dynamics
The aim of this tutorial is to help you getting started using the CPMD program. CPMD is an ab initio electronic structure and molecular dynamics (MD) program using a plane wave/pseudopotential implementation of density functional theory (DFT). It is mainly targeted at Car-Parrinello MD simulations, but also supports geometry optimisations, Born-Oppenheimer MD, path integral MD, response functions, excited states and calculations of some electronic properties. For further information you may want to take a look at the CPMD consortium homepage at www.cpmd.org.
Born-Oppenheimer (BO) MD: at each step of atomic MD, one solves the DFT SCF problem for a fixed atomic configuration, computes the wavefunction, and from it the (quantum mechanical) forces acting on the nuclei. This allows it to propagate the atoms classically to the next configuration. (In simple words, one computes energy and forces, move the atoms, re-compute energy and forces for the new configuration, move the atoms, and so on.)
Car-Parrinello: One propagates simultaneously a set of atoms and a fictitious set of electrons. From the fictitious set of electrons, one can estimate the forces on the atoms and move the atoms. Atoms and fictitious electrons have to be independent, and their dynamics is described using a generalised Lagrangian formalism, rather than simply propagating the atoms using Newton's equations like in Born-Oppenheimer MD. In the CPMD manual, it is written:
“The basic idea of the Car-Parrinello approach can be viewed to exploit the time-scale separation of fast electronic and slow nuclear motion by transforming that into classical-mechanical adiabatic energy-scale separation in the framework of dynamical systems theory. In order to achieve this goal the two-component quantum / classical problem is mapped onto a two-component purely classical problem with two separate energy scales at the expense of loosing the explicit time-dependence of the quantum subsystem dynamics.”
Born-Oppenheimer MD is easier to implement in practice, but slower. Also it requires a very well converged calculation of the forces at each step of dynamics. The Car-Parrinello only requires a very accurate ground state at the first step. The rest of the dynamics sustains itself, provided it is carried out in the right way.
If you would like to know more about the theoretical background of CPMD or BOMD, please read CPMD manual Sections 6.6 and 6.7.
2. Code Running Environment
2.1 Before running the job
Before you run CPMD it is a good idea to copy all you need in your local directory.
- Copy the folder lab from the usb disk to your laptop.
- open the terminal application
- open two windows, one will be for your local mac and the other for the hpc cluster
- on one of them login to the hpc cluster by typing ssh yourname@login.cx1.hpc.ic.ac.uk
- the cluster has two key directories /WORK/yourname and /HOME/yourname do all the calculations from the WORK/yourname directory, from now on this directory will be referred to simply as work
- cd to your work directory and make a new directory called liquids, cd into the new directory, this is where you will run all your jobs from
- copy the file from your desktop to hpc:
scp filename yourname@login.cx1.hpc.ic.ac.uk:/work/yourname/foldername/.
- unzip the tar file
tar -xvf filename.tar
2.2 How to run the job
In the directory cpmd you should have two subdirectories: clw1 and clw63.
To run CPMD type qsub runcpmd_wf (or other run scripts depending on the type of calculation such as runcpmd_md) in the clw1 or clw63 directory.
2.3 Using VMD to analyse results
For visualisation of the results you may want to take a look at the tutorial on visualising results from CPMD with the VMD program under www.theochem.ruhr-uni-bochum.de/go/cpmd-vmd.html
3. CPMD
3.1 Theoretical background
3.1.1 Pseudopotentials (read CPMD manual section 4.2)
3.1.2 Kohn-Sham orbitals, Wannier orbitals and Effective molecular orbitals
The EMO method for analysing the electronic structure of solutes and solvents starts by considering the Kohn-Sham Hamiltonian and energy eigenvalues obtained from the ab initio molecular dynamics. A band structure similar to that obtained for solid state systems emerges. While the Kohn-Sham Hamiltonian is diagonal, and the Kohn-Sham orbital energies are uniquely defined, the Kohn-Sham orbitals are delocalised Bloch states which extend over all space and the molecular picture is lost, as shown in Figure 1 (a). This makes it very difficult to analyse interactions between an ion and the solvating water molecules. Maximally localized Wannier orbitals are local to atomic centers. However, the Wannier Hamiltonian is not diagonal and thus there is no unique energy eigenvalue that can be associated with each orbital, as shown in Figure 1 (c). An intermediate representation where orbitals are local to a molecule, but are delocalised within the molecule has been developed. The EMOs are obtained by assigning Wannier orbitals to a specific molecule, and block-diagonalising the Wannier Hamiltonian within the molecular subspace, as shown in Figure 1 (b). This analysis has been successfully applied to pure liquid water. An isolated water molecule has C2v symmetry and valence MOs 1a1, 1b2, 2a1 and 1b1. A liquid environment is disordered and there is no symmetry. However, Figure 1 (b) shows that the EMOs strongly resemble the standard MOs for a single water molecule in vacuum. This close resemblance is by no means predetermined and is a validation of the EMO method. For clarity, we will continue to use the isolated atomic or molecular symmetry labels when we refer to the individual atomic orbitals, MOs and bands in the aqueous system.
3.2 Functional
We will use a gradient corrected functional (BLYP) instead of LDA throughout the tutorial as this combination of generalized gradient approximations to exchange and correlation has been shown to give good results for the structure and dynamics of water.
Figure 1 Wavefunction diagrams from a CPMD simulation of pure liquid water related to graphical representations of the respective Hamiltonian matrices: (a) Kohn-Sham Hamiltonian, (b) EMO Hamiltonian and (c) Wannier Hamiltonian.
3.3 Generating input files and understanding output files
The first example will demonstrate some of the basic steps of performing a CPMD calculation with a very simple system: 1 water molecule+Cl-, and a very simple task: calculate the electronic structure. We will use that as an example to have a look at the input file format, and how to read the output.
- emacs filename
- and then look inside and follow the instructions below
3.3.1 Wavefunction optimisation
(a) Input File format
For nearly all CPMD calculations, you first have to optimise the wavefunction of your system, and use that as a base for further calculations. For our first calculation you’ll need the input file inp.clw_blyp_wf and the pseudopotential files H_MT_BLYP.psp, O_MT_BLYP.psp, Cl_MT_BLYP.psp.
Now let’s have a look at the input file. The input file is organised in sections which start with &NAME and end with &END. Everything outside those sections is ignored. All keywords have to be in upper case or else they will be ignored. The sequence of the sections does not matter, nor does the order of keywords. A minimal input file must have a &CPMD, &SYSTEM and an &ATOMS section.
&CPMD OPTIMIZE WAVEFUNCTION CONVERGENCE ORBITALS 1.0d-7 &END
This first part of the &CPMD section instructs the program to do a wavefunction optimization (i.e. a single point calculation) with a tight convergence criterion (the default is 1.0d-5).
Exercise: look up in the CPMD manual about the usage of the keywords ODIIS, STORE, TIMESTEP, EMASS, ISOLATED MOLECULE.
&DFT FUNCTIONAL BLYP GC-CUTOFF 0.1D-06 &END
The &DFT section is used to select the density functional and related parameters. In this case we go with the BLYP functional.
&SYSTEM ANGSTROM SYMMETRY 0 CELL 11.800 1.0 1.0 0 0 0 CUTOFF 70.0 CHARGE -1.0 &END
The &SYSTEM section contains various parameters related to the simulation cell and the representation of the electronic structure. The keywords SYMMETRY, CELL and CUTOFF are required and define the (periodic) symmetry, shape, and size of the simulation cell, as well as the plane wave cutoff (i.e. the size of the basis set). The keyword Angstrom additionally indicates that all lengths and coordinates are given in angstrom (not in a.u.).
Exercise: look up in the CPMD manual about the usage of the keyword CHARGE.
Finally the &ATOMS section is needed to specify the atom coordinates and the pseudopotentials, that are used to represent them. The detailed syntax of the pseudopotential specification is a bit complicated and will not be needed nor discussed here. If you want to know more, please have a look at the Further details of the Input section of the CPMD manual (P222).
- To run, type: qsub runcpmd_wf,
The code will generate the following output files once the calculation is completed after a few minutes: log.clw_blyp_wf, GOMETRY., RESTART.. The purpose of the wavefunction optimisation is to generate the RESTART. file which is needed for further calculations.
(b) Output File Format
To start the calculation, which should be completed in less than a minute once the job starts running. The main output of the CPMD program is now in the file log.clw_blyp_wf. Let’s have a closer look at the contents of this file.
PROGRAM CPMD STARTED AT: Thu Jan 20 16:23:24 2011
SETCNST| USING: CODATA 2006 UNITS
THE INPUT FILE IS: /work/lge/wannier/test/inp.clw_blyp_wf
THIS JOB RUNS ON: cx1-50-2-1.cx1.hpc.ic.ac.uk
THE CURRENT DIRECTORY IS:
/tmp/pbs.5020787.cx1
THE TEMPORARY DIRECTORY IS:
/tmp/pbs.5020787.cx1
THE PROCESS ID IS: 13680
Here we have some technical information about the environment, where this job was run.
SINGLE POINT DENSITY OPTIMIZATION
PATH TO THE RESTART FILES: ./
GRAM-SCHMIDT ORTHOGONALIZATION
MAXIMUM NUMBER OF STEPS: 1000 STEPS
MAXIMUM NUMBER OF ITERATIONS FOR SC: 1000 STEPS
PRINT INTERMEDIATE RESULTS EVERY 10001 STEPS
STORE INTERMEDIATE RESULTS EVERY 10 STEPS
NUMBER OF DISTINCT RESTART FILES: 1
TEMPERATURE IS CALCULATED ASSUMING AN ISOLATED MOLECULE
FICTITIOUS ELECTRON MASS: 500.0000
TIME STEP FOR ELECTRONS: 5.0000
TIME STEP FOR IONS: 5.0000
CONVERGENCE CRITERIA FOR WAVEFUNCTION OPTIMIZATION: 1.0000E-07
WAVEFUNCTION OPTIMIZATION BY PRECONDITIONED DIIS
THRESHOLD FOR THE WF-HESSIAN IS 0.5000
MAXIMUM NUMBER OF VECTORS RETAINED FOR DIIS: 5
STEPS UNTIL DIIS RESET ON POOR PROGRESS: 5
FULL ELECTRONIC GRADIENT IS USED
SPLINE INTERPOLATION IN G-SPACE FOR PSEUDOPOTENTIAL FUNCTIONS
NUMBER OF SPLINE POINTS: 5000
This section now gives you a summary of the parameters read in from the &CPMD section, or their respective default settings.
EXCHANGE CORRELATION FUNCTIONALS
LDA EXCHANGE: SLATER (ALPHA = 0.66667)
LDA CORRELATION: LEE, YANG & PARR
[C.L. LEE, W. YANG, AND R.G. PARR, PRB 37 785 (1988)]
GRADIENT CORRECTED FUNCTIONAL
DENSITY THRESHOLD: 1.00000E-07
EXCHANGE ENERGY
[A.D. BECKE, PHYS. REV. A 38, 3098 (1988)]
PARAMETER BETA: 0.004200
CORRELATION ENERGY
[LYP: C.L. LEE ET AL. PHYS. REV. B 37, 785 (1988)]
*** DETSP| SIZE OF THE PROGRAM IS 8748/ 110552 kBYTES ***
>>>>>>>> CENTER OF MASS HAS BEEN MOVED TO CENTER OF BOX <<<<<<<<
***************************** ATOMS ****************************
NR TYPE X(bohr) Y(bohr) Z(bohr) MBL
1 Cl 11.199323 11.186938 9.180419 3
2 O 10.970667 11.067318 15.135702 3
3 H 10.558139 11.108892 13.264684 3
4 H 12.820897 11.171631 15.013437 3
****************************************************************
NUMBER OF STATES: 8
NUMBER OF ELECTRONS: 16.00000
CHARGE: -1.00000
ELECTRON TEMPERATURE(KELVIN): 0.00000
OCCUPATION
2.0 2.0 2.0 2.0 2.0 2.0 2.0 2.0
This part of the output tells you which and how many atoms and electrons are used, what functional and what pseudopotentials were used, and what the values of some related parameters are.
POISSON EQUATION SOLVER : HOCKNEY COULOMB SMOOTHING RADIUS : 1.593 SYMMETRY: SIMPLE CUBIC LATTICE CONSTANT(a.u.): 22.29877 CELL DIMENSION: 22.2988 1.0000 1.0000 0.0000 0.0000 0.0000 VOLUME(OMEGA IN BOHR^3): 11087.72952 LATTICE VECTOR A1(BOHR): 22.2988 0.0000 0.0000 LATTICE VECTOR A2(BOHR): 0.0000 22.2988 0.0000 LATTICE VECTOR A3(BOHR): 0.0000 0.0000 22.2988 RECIP. LAT. VEC. B1(2Pi/BOHR): 0.0448 0.0000 0.0000 RECIP. LAT. VEC. B2(2Pi/BOHR): 0.0000 0.0448 0.0000 RECIP. LAT. VEC. B3(2Pi/BOHR): 0.0000 0.0000 0.0448 REAL SPACE MESH: 120 120 120 WAVEFUNCTION CUTOFF(RYDBERG): 70.00000 DENSITY CUTOFF(RYDBERG): (DUAL= 4.00) 280.00000 NUMBER OF PLANE WAVES FOR WAVEFUNCTION CUTOFF: 54804 NUMBER OF PLANE WAVES FOR DENSITY CUTOFF: 438632
This part of the output presents the settings read in from the &SYSTEM section of the input file and some derived parameters. […]
(K+E1+L+N+X) TOTAL ENERGY = -31.48398404 A.U.
(K) KINETIC ENERGY = 19.91015701 A.U.
(E1=A-S+R) ELECTROSTATIC ENERGY = -26.99115337 A.U.
(S) ESELF = 28.92331533 A.U.
(R) ESR = 0.74536593 A.U.
(L) LOCAL PSEUDOPOTENTIAL ENERGY = -19.20863875 A.U.
(N) N-L PSEUDOPOTENTIAL ENERGY = 2.63857242 A.U.
(X) EXCHANGE-CORRELATION ENERGY = -7.83292134 A.U.
GRADIENT CORRECTION ENERGY = -0.43675318 A.U.
After some output to report the setup of the initial guess for the electronic structure, we now see a summary of the various energy contribution of the total energy of the system, based on the initial guess. Now the program is ready to start the wavefunction optimisation.
Starting from the initial guess based on atomic wavefunctions the wavefunction for the total system is now calculated with an optimisation procedure. You can follow the progress of the optimisation in the output file.
NFI GEMAX CNORM ETOT DETOT TCPU 1 4.556E-02 4.038E-03 -31.483984 0.000E+00 5.88 2 1.580E-02 1.381E-03 -32.031355 -5.474E-01 5.90 3 1.347E-02 6.554E-04 -32.139961 -1.086E-01 5.93 4 7.667E-03 3.432E-04 -32.162687 -2.273E-02 5.95 5 3.490E-03 1.256E-04 -32.168927 -6.240E-03 5.97 6 2.655E-03 6.618E-05 -32.170107 -1.180E-03 5.98 7 1.888E-03 3.959E-05 -32.170488 -3.805E-04 5.99 8 1.004E-03 2.399E-05 -32.170647 -1.597E-04 5.97 9 6.207E-04 1.431E-05 -32.170689 -4.167E-05 5.96 10 4.309E-04 8.966E-06 -32.170703 -1.373E-05 5.96
NFI: Step number (number of finite iterations) GEMAX: largest off-diagonal component CNORM: average of the off-diagonal components ETOT: total energy DETOT: change in total energy to the previous step TCPU: (CPU) time for this step
And you can see that the calculation stops after the convergence of 1.0d-7 has been reached for the GEMAX value.
****************************************************************
* *
* FINAL RESULTS *
* *
****************************************************************
****************************************************************
* ATOMIC COORDINATES *
****************************************************************
1 Cl 11.199323 11.186938 9.180419
2 O 10.970667 11.067318 15.135702
3 H 10.558139 11.108892 13.264684
4 H 12.820897 11.171631 15.013437
****************************************************************
ELECTRONIC GRADIENT:
MAX. COMPONENT = 9.69721E-08 NORM = 8.07825E-10
TOTAL INTEGRATED ELECTRONIC DENSITY
IN G-SPACE = 16.000000
IN R-SPACE = 16.000000
(K+E1+L+N+X) TOTAL ENERGY = -32.17071373 A.U.
(K) KINETIC ENERGY = 18.31556847 A.U.
(E1=A-S+R) ELECTROSTATIC ENERGY = -26.97530784 A.U.
(S) ESELF = 28.92331533 A.U.
(R) ESR = 0.74536593 A.U.
(L) LOCAL PSEUDOPOTENTIAL ENERGY = -18.46114493 A.U.
(N) N-L PSEUDOPOTENTIAL ENERGY = 2.17769037 A.U.
(X) EXCHANGE-CORRELATION ENERGY = -7.22751980 A.U.
GRADIENT CORRECTION ENERGY = -0.41653735 A.U.
****************************************************************
Here we have the final summary of the results from our single point calculation. Please note that regardless of the input units, coordinates in the CPMD output are always in atomic units.
Other Output files:
Apart from the console output, our CPMD run created a few other files. Most importantly the restart file RESTART., which contains the final state of the system when the program terminated. This is needed to start other calculations, which need a converged wavefunction as a starting point. The file GEOMETRY. contains the coordinates of the atoms in atomic unit. This can be converted to a .xyz file.
3.3.2 Car-Parrinello Molecular dynamics
Based on the previously calculated electronic structure, we can now start a Car-Parrinello Molecular dynamics calculation. Note that although you can start a CP-MD run from a non-converged wavefunction (e.g. by not restarting from a pre-optimised wavefunction), you will be far away from the Born-Oppenheimer surface, and thus your result will be unphysical.
To run, type:
cp RESTART. RESTART qsub runcpmd_md
3.3.2.1 Input for CP dynamics
For the CP-MD job you need a new input file, inp.clw_blyp_md, which should be copied into the same directory, where you started the wavefunction optimisation run. If you compare it to the previous input files, you will find, that the only changes are again only in the &CPMD section of the input files.
&CPMD MOLECULAR DYNAMICS CP RESTART WAVEFUNCTION COORDINATES VELOCITIES NOSEP LATEST QUENCH BO MIRROR CONVERGENCE ORBITAL 1.0D-7 ODIIS 5 MAXSTEP 3000 STORE 20 TIMESTEP 5 EMASS 500.0 NOSE IONS 298.0 1100 TRAJECTORY SAMPLE 10 ISOLATED MOLECULE &END
The keyword MOLECULAR DYNAMICS CP defines the job type. Furthermore we tell the CPMD program to pick up the previously calculated wavefunction and coordinates from the latest restart file (which is named RESTART by default). MAXSTEP limits the MD to 3000 steps and the equations of motion will be solved for a time step of 5 atomic units (~0.12 femtoseconds). The temperature of the system will be initialised to 298K via the NOSE IONS keyword with Nose thermostat.
- Now start the CPMD program once more: qsub runcpmd_md
This run should be completed in a few minutes.
3.3.2.2 CP dynamics output
The output of the CPMD program is now in the file log.clw_blyp_md. There are also some new files (TRAJECTORY, ENERGIES). We will have a look at the output file first.
CAR-PARRINELLO MOLECULAR DYNAMICS
PATH TO THE RESTART FILES: ./
RESTART WITH OLD ORBITALS
RESTART WITH OLD ION POSITIONS
RESTART WITH OLD VELOCITIES
RESTART WITH OLD ION THERMOSTAT
RESTART WITH LATEST RESTART FILE
ITERATIVE ORTHOGONALIZATION
MAXIT: 30
EPS: 1.00E-06
MAXIMUM NUMBER OF STEPS: 3000 STEPS
The header is unchanged up to the point where the settings from the &CPMD section are printed. As you can see, the program has recognised the RESTART and the MAXSTEP keywords. (NOTE: in the CPMD code atoms are sometimes referred to as ions. This is due to the pseudopotential approach, where you integrate the core electrons into the (pseudo)atom which then could also be described as an ion.)
TIME STEP FOR ELECTRONS: 5.0000
TIME STEP FOR IONS: 5.0000
QUENCH SYSTEM TO THE BORN-OPPENHEIMER SURFACE
TRAJECTORIES ARE SAVED ON FILE EVERY 10 STEPS
ELECTRON DYNAMICS: THE TEMPERATURE IS NOT CONTROLLED
ION DYNAMICS: TEMPERATURE CONTROL (NOSE-HOOVER THERMOSTATS)
TARGET TEMPERATURE(KELVIN): 2.980000E+02
CHARACTERISTIC FREQUENCY(CM**-1): 1100.00
NOSE PARAMETERS
NUMBER OF THERMOSTATS (IONS) : 3
NUMBER OF THERMOSTATS (ELECTRONS) : 3
NUMBER OF THERMOSTATS (CELL) : 3
SCALING FOR ELEC. DOF : 6.00
NUMBER OF YOSHIDA-SUZUKI STEPS : 7
NUMBER OF INTEGRATION CYCLES (NIT): 1
This part of the output tells us, that the TIMESTEP 5.0 keyword was recognised (which is the default timestep), that the trajectory will be recorded and that NOSE-HOOVER thermostat is used.
RESTART INFORMATION READ ON FILE ./RESTART
Here we get notified, that the program has read the requested data from the restart file.
NFI EKINC TEMPP EKS ECLASSIC EHAM DIS TCPU
1 0.00000 0.0 -32.17071 -32.17010 -32.17010 0.458E-10 3.25
2 0.00000 0.0 -32.17071 -32.17010 -32.17010 0.727E-09 3.24
3 0.00000 0.0 -32.17071 -32.17010 -32.17010 0.365E-08 3.24
4 0.00000 0.0 -32.17071 -32.17010 -32.17010 0.114E-07 3.24
5 0.00000 0.1 -32.17071 -32.17010 -32.17010 0.275E-07 3.24
6 0.00000 0.1 -32.17071 -32.17010 -32.17010 0.565E-07 3.25
7 0.00000 0.1 -32.17071 -32.17010 -32.17010 0.104E-06 3.24
8 0.00000 0.2 -32.17072 -32.17010 -32.17010 0.175E-06 3.24
9 0.00000 0.2 -32.17072 -32.17010 -32.17010 0.279E-06 3.25
10 0.00000 0.3 -32.17072 -32.17010 -32.17010 0.422E-06 3.25
11 0.00000 0.3 -32.17072 -32.17010 -32.17010 0.614E-06 3.25
12 0.00000 0.4 -32.17072 -32.17010 -32.17010 0.865E-06 3.25
13 0.00000 0.4 -32.17072 -32.17010 -32.17010 0.119E-05 3.24
14 0.00000 0.5 -32.17072 -32.17010 -32.17010 0.159E-05 3.24
15 0.00000 0.6 -32.17072 -32.17010 -32.17010 0.209E-05 3.25
16 0.00000 0.6 -32.17072 -32.17010 -32.17010 0.269E-05 3.25
17 0.00000 0.7 -32.17072 -32.17010 -32.17010 0.343E-05 3.24
18 0.00000 0.8 -32.17072 -32.17010 -32.17010 0.429E-05 3.25
19 0.00000 0.9 -32.17072 -32.17010 -32.17010 0.532E-05 3.24
20 0.00000 1.0 -32.17072 -32.17010 -32.17010 0.652E-05 3.24
...
2981 0.00015 274.5 -32.16644 -32.17025 -32.17010 0.443E+01 3.25
2982 0.00015 268.3 -32.16638 -32.17025 -32.17010 0.444E+01 3.25
2983 0.00014 262.2 -32.16632 -32.17025 -32.17010 0.444E+01 3.26
2984 0.00014 256.0 -32.16627 -32.17024 -32.17010 0.445E+01 3.25
2985 0.00014 249.9 -32.16621 -32.17024 -32.17010 0.446E+01 3.26
2986 0.00013 243.8 -32.16615 -32.17024 -32.17010 0.446E+01 3.25
2987 0.00013 237.7 -32.16610 -32.17023 -32.17010 0.447E+01 3.25
2988 0.00013 231.7 -32.16604 -32.17023 -32.17010 0.447E+01 3.25
...
2997 0.00010 179.1 -32.16554 -32.17020 -32.17010 0.452E+01 3.25
2998 0.00010 173.6 -32.16549 -32.17020 -32.17010 0.452E+01 3.26
2999 0.00009 168.1 -32.16544 -32.17020 -32.17010 0.453E+01 3.26
3000 0.00009 162.6 -32.16539 -32.17019 -32.17010 0.453E+01 3.25
After some more output, we already discussed for the wavefunction optimisation, this is now part of the energy summary for a Car-Parrinello-MD run.
NFI: Step number (number of finite iterations) EKINC: (fictitious) kinetic energy of the electron (sub-)system TEMPP: Temerature (=kinetic energy / degrees of freedom) for atoms (ions) EKS: Kohn-Sham Energy: equivalent to the potential energy in classical MD ECLASSIC: Equivalent to the total energy in a classical MD (ECLASSIC = EHAM-EKINC) EHAM: total energy, should be conserved DIS: mean squared displacement of the atoms from the initial coordinates TCPU: (CPU) time needed for this step
To read about why do the MD part and how to analyse results, please read this CPMD paper on water.
Reference: Kuo et al, J. Phys. Chem. B 2004, 108, 12990-12998
Note: For a meaningful Car-Parrinello MD, EKINC has to be (and stay) very small (although for larger systems with more electrons, the absolute value of EKINC will be larger, i.e. no drift in EKINC.
Exercise: plot EKINC vs NFI, TEMPP vs NFI, EKS vs NFI, ECLASSIC vs NFI, EHAM vs NFI. Discuss the trend of these curves.
MEAN VALUE +/- RMS DEVIATION
<x> [<x^2>-<x>^2]**(1/2)
ELECTRON KINETIC ENERGY 0.000125 0.105171E-03
IONIC TEMPERATURE 233.0551 194.923
DENSITY FUNCTIONAL ENERGY -32.169126 0.152974E-02
CLASSICAL ENERGY -32.170228 0.105180E-03
CONSERVED ENERGY -32.170104 0.247772E-05
NOSE ENERGY ELECTRONS 0.000000 0.00000
NOSE ENERGY IONS -0.003317 0.200165E-02
CONSTRAINTS ENERGY 0.000000 0.00000
RESTRAINTS ENERGY 0.000000 0.00000
ION DISPLACEMENT 1.25973 1.33130
CPU TIME 3.2487
Finally we get a summary of some averages and root mean squared deviations for some of the monitored quantities. This is quite useful to detect unwanted energy drifts or too large fluctuations in the simulation.
Exercise: visualise the motion of the system of (1 water + Cl-)
Type: ./traj2xyz.x < TRAJECTORY.clw_blyp_md > TRAJ.xyz
You can load the file TRAJ.xyz directly into the molecular visualisation program VMD.
Launch VMD program File -> Load files for: New Molecule click ‘Browse’, choose the file TRAJ.xyz click ‘Load’
The trajectory file can be used to generate structural properties such as radial distribution functions and analyse the solvation shells such as Cl-O distances.
3.3.3 Kohn-Sham energies
Altogether the first part of the CPMD input (file name: inp.clw_blyp_ksl) now contains:
&CPMD RESTART WAVEFUNCTION VELOCITIES COORDINATES LATEST KOHN-SHAM ENERGIES 0 MAXSTEP 1000 STORE 10 TIMESTEP 5.0 EMASS 500.0 ISOLATED MOLECULE &END
Output file (log.clw_blyp_ksl)
EIGENVALUES(EV) AND OCCUPATION:
1 -19.9458449 2.00000000 2 -11.8998056 2.00000000
3 -7.3120943 2.00000000 4 -4.5277656 2.00000000
5 -2.2061693 2.00000000 6 -0.5375340 2.00000000
7 -0.3936740 2.00000000 8 -0.3633460 2.00000000
CHEMICAL POTENTIAL = -0.3633465899 EV
This part lists the Kohn-Sham energies.
3.3.4 Kohn-Sham Orbitals
Altogether the first part of the CPMD input (file name: inp.clw_blyp_ksorb) now contains:
&CPMD KOHN-SHAM ENERGIES 8 RESTART WAVEFUNCTION COORDINATES LATEST DAVIDSON DIAGNOALIZATION RHOOUT BANDS 8 -1 -2 -3 -4 -5 -6 -7 -8 &END
This generages a series of files with the names WAVEFUNCTION.1, WAVEFUNCTION.2, …, WAVEFUNCTION.8 which have to be converted to cube format with cpmd2cube.x -o clw1 WAVEFUNCTION.1
3.3.5 Wannier Analysis
To calculate the Wannier centers or orbitals for our system we need to do a properties calculation starting from the previously generated wavefunction.
For the Wannier orbitals the keywords LOCALIZE and WANNIER WFNOUT are required. This generates a series of files with the names WANNIER_1.* which have to be converted to cube format with cpmd2cube.x -o clw1 WANNIER_1.1.
Altogether the first and second parts of the CPMD input to calculate Wannier enters now contains:
&CPMD RESTART WAVEFUNCTION LATEST PROPERTIES WANNIER DOS &END &PROP LOCALIZE &END
To plot all the WANNIER ORBITALS, replace the first part with:
&CPMD RESTART WAVEFUNCTION COORDINATES LATEST PROPERTIES WANNIER WFNOUT ALL &END
Exercise: use VMD to visualise the Kohn-Sham and Wannier orbitals. Describe the Kohn-Sham and Wannier Orbitals of the system. (Hint: One of the Kohn-Sham orbitals looks like this. Discuss delocalisation/localisation.)
Figure 2 Kohn-Sham Orbital
Note: Different phases can be visualised in VMD by creating two representations from the same data set and just using isovalues with opposite sign for each of them. The input and cube files used in this section are:
WAVEFUNCTION.* or WANNIER_1.*
clw1.*.cube (for Kohn-Sham orbitals) or clw11.*.cube (for Wannier orbitals)
VMD:
-> Launch VMD program -> File -> Load files for: New Molecule -> click ‘Browse’, choose the file TRAJ.xyz -> click ‘Load’ -> Graphics -> Representations -> Drawing Method -> Choose CPK -> repeat the above procedure and load the cube file again but this time choose Isosurface as Drawing method and play around with Isovalues.
To change colour of the background:
-> Graphics -> Colors -> Display -> Background -> Choose white
4. Exercises (63 water + Cl-)
The next step is a more typical (while ambitious) application of the CPMD code: a Car-Parrinello MD simulation of a bulk system with water and a chloride ion. In this specific example, we try to look at how the electronic structure changes when a chloride ion is solvated in liquid water. Our system will consist of 63 water molecules and one chloride ion (note the CHARGE keyword in the &SYSTEM section). To speed up the equilibration phase, we start from a restart configuration that has been equilibrated with classical MD for about 3 ns using the SPC/E water potential and CPMD for 31.56 ps. Since water has a dipole moment, you have to keep in mind, that we are calculating a system with periodic boundary conditions, so a water molecules ‘sees’ its images and interacts with them. There are methods implemented in CPMD to compensate this effect, but we won’t discuss them here.
Now run jobs to get data.
4.1 Run a wavefunction optimisation (use inp.w63Cl-wf).
4.2 We want to run the MD at 298 Kelvin, so now run a short MD with temperature rescaling for the atoms and no thermostat for the electrons (use inp.w63Cl-md)
4.3 Now you have your MD trajectory files, analyse them
4.4 Liquid water structure and hydrogen bonding
4.4.1 Generate a movie of trajectories using VMD which will read in the TRAJ.xyz file.
4.4.2 what to look at => Radial Distribution Functions
4.4.3 Compare your traj to the one we have made
4.4.4 Do you see water molecules entering and exiting solvation shell?
4.5 Plotting Kohn-Sham energies, Kohn-Sham orbitals, Wannier centers and Wannier orbitals
4.6 Follow the instructions above for clw1 but with your clw63 file (do it for last step of your trajectory) in reality we would want to examine 50-100 points sampled from the trajectory
4.7 Each student to look at one of the following:
4.7.1 Plot density of states for KS, KS orbitals
4.7.2 Plot Wannier orbitals and centers
4.7.3 Plot Cl-water distances over time