Talk:Mod:Hunt Research Group/wannier centre
Characterising the electronic states of Cl- in 63 water molecules (research notes by Ling)
1. Wannier Centres
My results of [63 water + Cl-] at 1.32 ps.
Fig. 1 Wannier centres of [63 water molecules + Cl-]. Blue - Wannier centers, Red - Oxygen, Grey - Hydrogen, Cyan - Chlorine ion.
Each oxygen atom is in the centre of the Wannier center pyramid. Each Wannier centre on a water molecule correspond to a sp3 hybrid orbital. Two are the O-H bonds, and two the lonepairs. Cl- is also in the centre of the Wannier center pyramid as there are 4 occupied valence orbitals on Cl- (3s and 3 3p). Evidently they hybridise too, to give four sp3 orbitals, in a tetrahedral arrangement.
2. KS energies and DOS
3. RDF
Appendix:
1. How to visualise Wannier Centres and atom coordinates
Method A (I recommend):
(1) Use the keyword 'WANNIER REFERENCE' followed by three parameters which are - half of cell size in bohr. E.g.
&CPMD PROPERTIES RESTART WAVEFUNCTION COORDINATES WANNIER REFERENCE -11.7352 -11.7352 -11.7352 &END &PROP LOCALIZE &END
&DFT NEWCODE FUNCTIONAL BLYP GC-CUTOFF 0.2D-06 &END
&SYSTEM SYMMETRY 1 CELL 23.4704 1.0 1.0 0 0 0 CUTOFF 70.0 CHARGE -1 &END
Note: During visualisation, make sure there are no bonds between Wannier centers and hydrogen atoms.
(2) Visualise (using VMD) the second set of coordinates in IONS+CENTERS.xyz, which include the coordinates of Wannier centers after localisation.
Method B:
I used coordinates of both atoms and Wannier centres from output, and put them back to the same cell.
(1) Use the coordinates of Wannier Centres after the optimisation in the output file (e.g. log.w63Cl-wanorb). Note: the coordinates in the output file are in bohr, need to convert to angstrom by 1 Bohr = 0.52917720859 Angstrom, using 'awk', e.g.
awk '{printf "%.12f %.12f %.12f \n", $1*0.5291772108, $2*0.5291772108, $3*0.5291772108}' w63Cl-wannier_center.xyz > 1.xyz
(2) Use the coordinates of atoms in the output file (e.g. log.w63Cl-wanorb). Note: the coordinates in the output file are in bohr, need to convert to angstrom by 1 Bohr = 0.52917720859 Angstrom. Put both sets of coordinates in a file called w63Cl-wannier_center_atoms.xyz
(3) Code to put Wannier Centres and atoms to the same cell (convert.f90). Use ./a.out < w63Cl-wannier_center_atoms.xyz, then I would get a file called output.xyz which can be visualised by VMD (please refer to my VMD notes on how to visualise and change the colour of atoms).
PROGRAM COORDINATES
implicit none integer, parameter :: NumAtoms=446 real, dimension(3) :: atoms character*2 :: symbols integer :: i, j real :: l character*30 :: ch
!cell size = 23.4704 Angstrom l = 23.4704
open(unit=7,file="output.xyz") read(5,*) ch write(7,*) ch read(5,*) ch write(7,*) ch
!The coordinates are in Angstrom
do i = 1, NumAtoms
read(5,*) symbols, atoms
atoms = atoms
if (trim(symbols)=="D") then
atoms = atoms + l/2.
end if
write(7,'(a,3f20.6)') symbols, atoms
end do
close(7)
END PROGRAM COORDINATES
2. How to find which KS energies are ion HOMO and which ones are water band?
First look at the density of states, and find the Fermi energy. Most likely the Cl- states will be just below the Fermi energy, in the water gap. Once we know the possible ion states (can visualise the corresponding orbitals), and see if they are localised on the ion.
Note: In ab inition MD we are using the Born-Oppenheimer approximation, so whether the atoms are thermodynamically in equilibrium or not does not make any difference for the electronic problem. Because the electron always see the nuclei as frozen in some fixed configuration.
It might well be that there are many empty states at zero energy. But take into account that the absolute values of the eigenvalues have no meaning, because the system is infinite. So all the eigenvalues need first to be referred to the Fermi energy.
If you compute a density of states, and then compare to some published results, you will see immediately if everything is OK.
My results at 1.32ps:
251 -2.4813154 2.00000000 252 -2.3225673 2.00000000 253 -2.2307503 2.00000000 254 -2.1896318 2.00000000 255 -2.1063431 2.00000000 256 -1.9922699 2.00000000 257 0.2844580 0.00000000 258 0.4295585 0.00000000 259 0.9961278 0.00000000 260 1.9684643 0.00000000 261 3.1466093 0.00000000 262 3.2006318 0.00000000 263 3.3180696 0.00000000 264 3.3661100 0.00000000 265 3.5236525 0.00000000 266 3.5883758 0.00000000 267 4.0884711 0.00000000 268 4.1854061 0.00000000 269 4.2129189 0.00000000 270 4.3560058 0.00000000 271 4.3683808 0.00000000 272 4.3995624 0.00000000 273 4.5195778 0.00000000 274 4.5428412 0.00000000 275 4.5786353 0.00000000 276 4.6169247 0.00000000 277 4.6728144 0.00000000 278 4.8055353 0.00000000 279 5.0279848 0.00000000 280 5.1143694 0.00000000 281 5.2572429 0.00000000 282 5.3395347 0.00000000 283 5.3693926 0.00000000 284 5.4557988 0.00000000 285 5.5166133 0.00000000 286 5.5705186 0.00000000 287 5.7131019 0.00000000 288 5.8420968 0.00000000 289 5.9221986 0.00000000 290 5.9497781 0.00000000 291 6.0387123 0.00000000 292 6.0945736 0.00000000 293 6.1667204 0.00000000 294 6.2048787 0.00000000 295 6.2851450 0.00000000 296 6.3167766 0.00000000 CHEMICAL POTENTIAL = -1.9922706859 EV
256 - highest occupied state. The chemical potential here is taken to be the energy of the HOMO. In general, the Fermi energy is the midway between the HOMO and the LUMO (i.e. at the centre of the gap).
The values make sense.
I have computed 40 empty states, to get a reasonable description of the conduction band of water as well.
To do: compute DOS. (Manual P182: KPERT, P185 LDOS, P219 WANNIER DOS).
Note: Manual P245
"Q: I have a problem with visualising unoccupied orbitals. When I use RHOOUT BANDS or CUBEFILE ORBITALS after the wavefunction optimization I get only occupied orbitals. If I add one empty state when optimizing wavefunction the program never reaches convergence.
A: The most efficient way to calculate unoccupied orbitals is to firrst optimize the occupied orbitals and then restart the calculation using the run option
KOHN-SHAM ENERGIES n
where n ist the number of unoccupied orbitals. This will diagonalize the Kohn-Sham Potential (defined by the occupied orbitals alone).
To test if everything goes fine, you can check the total energy printed at the beginning of this job, it should be exactly the one at the end of the optimization. In addition, if you don't change the default convergence criteria, the number of converged Kohn-Sham states should be equal to the number of occupied states in the first step."