In this exercise, we will look at a more complicated example than in Exercise 1, namely liquid water. We start by computing the PDOS and compare it the single H$_2$O molecules. This time we will have many more orbitals so in order to visualize our results, we will apply a Gaussian broadening to the results. After, we will do a short Born-Oppenheimer MD run and learn how the basic input and output.
In this part of the exerise we will compute the PDOS. Our model water will contain 64 H$_2$O molecules in a 12 Å box under periodic boundary conditions. The input is very similar to the first exercise, but since we now have a much more coordinates, it is convinient to import them from a separate file. Below is a an example input you can use to perform a single-point calculation extracting the PDOS.
&GLOBAL
PROJECT liquid_water
RUN_TYPE ENERGY
IOLEVEL LOW
&END GLOBAL
&FORCE_EVAL
METHOD QS
&DFT
BASIS_SET_FILE_NAME GTH_BASIS_SETS
POTENTIAL_FILE_NAME GTH_POTENTIALS
LSD 0
&QS
METHOD GPW
EPS_DEFAULT 1.0E-10
&END QS
&SCF
MAX_SCF 300
EPS_SCF 1.0E-06
SCF_GUESS ATOMIC
&MIXING
METHOD DIRECT_P_MIXING
ALPHA 0.45
&END MIXING
&DIAGONALIZATION
ALGORITHM STANDARD
&END DIAGONALIZATION
&PRINT
&RESTART
BACKUP_COPIES 0
&END
&END PRINT
&END SCF
&MGRID
NGRIDS 4
CUTOFF 300
REL_CUTOFF 60
&END
&XC
&XC_FUNCTIONAL BLYP
&END XC_FUNCTIONAL
&END XC
&PRINT
&PDOS
FILENAME ./liquid_water
&END PDOS
&END PRINT
&END DFT
&SUBSYS
&CELL
ABC 12.4170 12.4170 12.4170
ANGLES 90.00 90.00 90.00
&END CELL
&TOPOLOGY
COORD_FILE_FORMAT XYZ
COORD_FILE_NAME liquid_water.xyz
CONNECTIVITY OFF
&END TOPOLOGY
&KIND H
BASIS_SET DZVP-GTH
POTENTIAL GTH-BLYP-q1
&END KIND
&KIND O
BASIS_SET DZVP-GTH
POTENTIAL GTH-BLYP-q6
&END KIND
&END SUBSYS
&END FORCE_EVAL
&MOTION
&GEO_OPT
OPTIMIZER BFGS
MAX_DR 3.0E-03
MAX_FORCE 4.5E-04
RMS_DR 1.5E-03
RMS_FORCE 3.0E-04
&END
&PRINT
&RESTART
BACKUP_COPIES 0
&END
&RESTART_HISTORY OFF
&END
&END PRINT
&END MOTION
To run this input, you also need to download the following coordinate file:
192
12.4170 12.4170 12.4170 90.0 90.0 90.0
O 0.636377 -1.361001 -3.178429 167.614871
O -6.713644 -0.661735 0.075884 16.831026
O 1.181465 8.391928 -6.312263 544.083738
O 6.225075 4.985516 -4.102950 223.150040
O 1.077589 0.157722 -0.792948 158.588916
O -3.192448 2.665559 1.269438 78.250523
O 4.186181 -4.648630 2.064529 363.121018
O 1.084065 -2.622383 1.720973 100.988895
O 3.961304 -3.500551 6.395341 207.317947
O -2.426954 -0.116530 -0.541931 80.001094
O 2.849800 3.072931 4.515085 246.847962
O -3.778861 -1.516854 3.044998 142.123443
O -1.780897 0.255074 -5.627996 108.140843
O -5.860061 1.510328 1.697173 230.500524
O -6.058608 6.103516 1.672394 341.329523
O -0.548833 8.593670 -3.964855 6.763006
O -6.787421 -5.431829 5.279242 69.756053
O 6.197776 3.761218 0.111023 314.505579
O -2.581713 2.191784 -2.246646 223.679673
O -0.358882 1.215490 1.611481 303.485510
O -3.894738 4.014102 -0.997240 289.269855
O 3.399572 -1.589453 -2.761965 683.549227
O -1.688429 -5.449113 3.709715 44.524651
O 6.161442 0.256797 -2.782453 534.772107
O -5.264801 -1.884111 5.118901 670.390824
O 2.519912 -5.437210 -0.035399 79.492377
O -0.071818 -0.834811 3.475668 264.670150
O 2.497556 2.421502 -4.060662 62.209680
O 1.779316 7.687595 3.186328 380.984079
O 0.201776 -7.017557 3.062847 67.106313
O 4.046782 -5.830823 -4.119788 69.282060
O -3.884227 -2.931983 0.875670 427.118732
O -3.202479 -0.955125 -3.113453 149.707977
O -4.764760 3.736065 -6.155987 570.261328
O -3.118051 -2.729306 -5.577398 587.438331
O 3.447757 -0.758854 1.726687 177.872630
O 4.446184 1.319425 5.829542 53.674527
O -1.037498 4.896488 -5.368932 411.249081
O 0.730083 0.435771 -5.007765 33.602971
O 2.151186 2.160608 1.918089 279.866499
O 5.343458 -0.919645 7.289043 609.167560
O 2.916396 4.299320 0.021636 127.674108
O 5.057850 4.583619 3.577332 163.917568
O -1.222523 3.286228 4.311124 27.659504
O -2.368688 2.725897 -4.812665 35.360124
O -2.897983 0.694401 4.307643 1095.553210
O -5.642040 -3.863037 -0.942542 82.182860
O -4.251798 6.263271 5.330620 104.803196
O -0.920362 -6.546462 0.438624 145.425559
O 4.804521 -2.394666 3.417408 221.285926
O -3.563471 -5.503069 1.323223 284.142419
O -1.141425 -2.876637 4.924695 388.984345
O 5.022231 2.714247 -2.863996 709.439925
O 7.172531 -2.749657 -3.169540 113.471231
O 1.334265 2.909853 -1.554424 295.885454
O -0.555670 4.921867 -1.944234 19.399452
O -2.900763 -5.736063 -4.566719 152.286243
O 2.075274 -0.389547 4.980292 85.557967
O 0.433354 -4.843016 -1.593907 972.475430
O -4.493104 -6.367184 -2.373044 47.395578
O -5.653472 1.651287 4.527714 46.538783
O 4.074561 -4.179597 -1.861619 84.531142
O -0.744505 -2.264736 -0.663989 304.470338
O 1.732339 4.835197 7.214335 73.739620
H 1.618291 -1.526563 -3.117950 420.352964
H 0.189990 -2.212772 -3.487410 421.321498
H -6.812885 -0.545112 -0.930279 1500.313802
H -6.251423 -1.525407 0.168724 304.046518
H 0.602180 9.155071 -6.524149 489.340001
H 2.133776 8.707665 -6.254781 154.380517
H 6.786086 5.400776 -3.364274 301.189381
H 6.789733 4.655344 -4.849357 90.904129
H 0.409178 0.306341 -0.110352 25.688680
H 0.639616 -0.355177 -1.494844 395.986768
H -4.127857 2.266753 1.501633 617.006760
H -2.987917 3.071299 2.148938 36.220243
H 4.277268 -3.896507 2.659175 290.035987
H 3.419083 -5.114309 2.528839 388.124547
H 0.509766 -1.911135 2.250576 1073.718253
H 1.878747 -2.066351 1.485318 805.710007
H 4.373583 -2.688703 6.749242 280.202224
H 4.711646 -4.085036 5.964371 428.989299
H -1.642436 -0.756828 -0.441072 485.253172
H -2.575023 0.404155 0.299777 201.028900
H 3.543736 3.634179 4.077058 858.394161
H 2.283813 3.716809 4.958762 753.660875
H -3.492489 -1.939537 3.866639 29.721930
H -3.469202 -0.579556 3.262107 79.927207
H -2.050900 -0.267623 -4.807888 238.525484
H -0.788270 0.260305 -5.525790 251.885934
H -6.174095 2.367261 1.381071 70.802318
H -5.779580 0.856643 0.974853 179.923934
H -6.632842 6.836133 1.868108 183.372371
H -5.121262 6.478161 1.616588 422.345445
H 0.192402 8.437457 -4.610564 68.949566
H -0.091114 8.306460 -3.114963 154.775180
H -6.956882 -5.101230 4.382713 38.346942
H -5.913658 -5.803664 5.187057 373.441453
H 6.245753 4.671478 0.474500 183.389450
H 7.153254 3.756662 -0.360072 429.366854
H -2.524827 1.245056 -1.959974 782.588879
H -2.470732 2.219782 -3.205491 10.546614
H 0.530770 1.689930 1.786587 53.974904
H -1.082239 1.845978 1.619624 184.996815
H -3.511880 3.280692 -1.605118 306.349695
H -3.451479 3.812062 -0.148594 351.607275
H 4.032795 -1.428993 -3.498578 225.184137
H 3.605797 -2.510895 -2.484469 560.187337
H -1.364239 -4.660963 4.135237 30.810044
H -2.545567 -5.578928 4.192120 210.247996
H 7.152896 0.331224 -2.764440 248.704160
H 5.852728 1.155669 -2.865559 168.479967
H -5.781592 -1.925170 4.276163 959.449833
H -5.903851 -1.584474 5.792820 667.400569
H 2.684682 -4.928776 0.743960 17.197726
H 2.769357 -6.378304 -0.003484 200.373806
H -0.418983 -0.067368 2.898211 140.563361
H 0.750496 -0.565436 4.038819 374.602009
H 2.069173 2.467111 -3.143679 624.302876
H 3.306855 1.980762 -3.904260 586.687014
H 1.364705 7.879493 4.109809 54.748087
H 1.341118 8.386086 2.584888 419.743650
H 0.928316 -6.358074 2.973841 707.660162
H -0.648802 -6.519197 3.322345 141.374271
H 4.847046 -6.429267 -4.168733 462.365583
H 4.284171 -5.184255 -4.810085 305.753847
H -3.154654 -2.790581 0.265587 452.650616
H -3.927884 -2.310131 1.639024 612.464607
H -4.018770 -1.618380 -3.076116 408.538263
H -2.859839 -0.864321 -2.194437 168.872836
H -5.191711 3.209094 -6.876638 433.874186
H -3.936247 3.321345 -5.869924 495.554644
H -3.868924 -2.553424 -6.149877 322.885381
H -3.206105 -2.187991 -4.819678 95.014812
H 3.020984 0.105586 1.759708 246.253249
H 4.048833 -0.684224 0.915677 134.674214
H 5.208767 1.365264 5.257056 83.884966
H 3.814664 2.035442 5.536889 552.642546
H -0.078327 4.949498 -5.157763 346.750021
H -1.597700 5.629327 -5.007949 959.620880
H 1.319507 1.194678 -4.854496 451.441577
H 0.820661 -0.249716 -4.293210 181.691634
H 2.325918 2.562747 2.815404 118.853508
H 2.410239 2.834726 1.148382 112.263781
H 4.954270 -0.104134 6.915226 81.834727
H 5.965311 -0.553463 7.979837 426.977783
H 2.442875 3.711271 -0.595321 3.045300
H 3.835014 4.110994 -0.237118 20.245645
H 5.223277 5.258621 2.869490 211.898078
H 5.092684 5.124820 4.408139 374.919057
H -1.448263 3.595627 5.210353 215.628060
H -0.416965 3.864827 4.037672 116.680107
H -1.702159 3.414418 -5.108024 26.187727
H -2.151180 1.799521 -5.229885 819.613090
H -2.563739 0.455339 5.185547 132.365647
H -2.519593 1.597205 4.245806 323.309625
H -6.570841 -4.114938 -0.625894 119.141393
H -5.249567 -3.348649 -0.194238 531.193434
H -4.006841 6.620008 6.170298 468.509495
H -4.385336 5.265871 5.552889 214.103176
H -1.664101 -6.065255 0.787609 501.768013
H -0.393868 -6.855340 1.267474 397.297206
H 4.355396 -2.424213 4.269310 833.851194
H 4.175551 -1.797993 2.878018 727.696780
H -2.939954 -5.418336 2.106708 61.831756
H -3.820096 -4.582157 1.173209 541.098029
H -1.746701 -2.440205 5.562455 147.181962
H -0.871309 -2.147110 4.339277 194.019066
H 5.113253 2.939116 -1.880597 49.617582
H 5.334314 3.531808 -3.408691 779.238718
H 7.055546 -3.496503 -3.797032 350.147677
H 6.923886 -3.109040 -2.312507 106.142905
H 0.960261 2.084917 -1.185341 63.912112
H 0.667000 3.653207 -1.686012 158.735414
H -1.170995 4.301608 -2.335690 56.804885
H -1.005102 5.018841 -1.073550 109.443123
H -2.273834 -5.003007 -4.425863 248.945606
H -3.367701 -5.839040 -3.748665 126.770794
H 1.600029 -0.105762 5.792579 51.368546
H 2.880005 0.144513 4.801890 223.542929
H 1.161929 -4.813603 -0.915777 101.363482
H 0.256593 -5.780733 -1.752487 444.325657
H -4.799807 -5.599160 -1.812579 250.199262
H -4.099071 -7.011965 -1.735563 454.266378
H -4.929740 0.962064 4.626477 363.561307
H -5.898601 1.585582 3.571661 177.267410
H 4.175887 -4.731434 -2.734608 132.314472
H 3.464890 -4.697521 -1.257727 318.937590
H -0.298830 -2.522941 0.145305 410.993119
H -0.664084 -3.102339 -1.226943 472.136378
H 2.190494 5.636486 7.493926 442.460984
H 2.308703 4.123715 7.631363 224.342520
The new section governing the reading of external coordinate files is
&TOPOLOGY COORD_FILE_FORMAT XYZ COORD_FILE_NAME liquid_water.xyz CONNECTIVITY OFF &END TOPOLOGY
where we specify the filename, format, and that we should not divide our atoms into molecules. When you have finished running the above input, try to use the following python script to convolve your spectra with Gaussian functions. You do this by first loading Python3 by typing
module load Anaconda3/5.0.1
and then typing
python pdos_conv.py
Note that you might have to change the parameters at the top of the file to get pleasing results. If a new window with the plot does not open, it is also saved in the current directory as “pdos.png”.
import numpy as np
import os
import matplotlib as mpl
mpl.use('Agg')
import matplotlib.pyplot as plt
# Parameters
FWHM = 2.5 # Full width half maximum
Emin=-16. # Starting energy
Emax= 4. # -----
DE = 0.01 # Fineness of the energies to convolve to
dE=np.arange(Emin,Emax,DE)
E_shift = 0
filename = 'liquid_water-k1-1.pdos'
def plot_pdos(txtfile):
if not os.path.isfile(txtfile):
print("Warning: File does not exist!")
return np.zeros(np.shape(dE)[0])
PDOS = np.asarray(open(txtfile).read().split())
if 'f' in PDOS:
start = 26
Norb = 4+3
elif 'd' in PDOS:
start = 25
Norb = 3+3
elif 'p' in PDOS:
start = 24
Norb = 2+3
else:
start = 23
Norb = 1+3
L = int(PDOS[-Norb])
pdos = PDOS[start:]
pdos = np.asarray(np.split(pdos,L))
E = pdos[:,1].astype(np.float)*27.2114 #a.u. -> eV
pdos = pdos[:,3:].astype(float)
s, p = convolve(E,pdos)
plt.plot(dE+E_shift,s, c='blue',linestyle='-',linewidth=2., label='O s')
plt.plot(dE+E_shift,p, c='red',linestyle='-',linewidth=2., label = 'O p')
plt.yticks([])
plt.ylabel('PDOS (arbitrary units)')
plt.xlabel('Binding energy (eV)')
plt.xticks(np.arange(-15,5,5),-np.arange(-15,5,5))
plt.xlim([-16,3])
plt.title('FWHM = ' + str(FWHM) + ' eV')
plt.legend(loc = 1)
plt.savefig('pdos.png')
plt.show()
def convolve(E,PDOS):
sigma = 1 / (2 * np.sqrt(2 * np.log(2) ) ) * FWHM
conv = np.zeros((len(dE),2))
for j in range(conv.shape[1]):
W = PDOS[:,j]
for i,e in enumerate(E):
conv[:,j] += np.exp(-(dE-e)**2 / (2 * sigma**2))*W[i]
return conv[:,0], conv[:,1]
plot_pdos(filename)
We will now attempt to do a short Born-Oppenheimer MD run on our system. The majority of the input from Part 1 can be reused while changing the run type to
RUN_TYPE MD
in the &GLOBAL section as well as adding the section
&MOTION
&MD
ENSEMBLE NVT
STEPS 15
TIMESTEP 1
TEMPERATURE 300.0
&THERMOSTAT
REGION MASSIVE
TYPE CSVR
&CSVR
TIMECON 5
&END CSVR
&END THERMOSTAT
&END MD
&END MOTION
The parameters in the MD section has been selected to finish in a few minutes, but feel free to experiment with different values and get a feel for how they work. To not print too much, also modify the &PDOS section to
&PDOS
FILENAME ./liquid_water
APPEND
&EACH
MD 3
&END
&END PDOS
which means that we only print the PDOS every third step, as well as appending the spectra instead of overwriting them.
When you have finished the calculation, you should have created some new files, including
Examine liquid_water-1.ener and make sure that the temperature and total energy remains reasonably stable during the run. As a final exercise, use a slightly modified python code from Part 1 to see how the PDOS spectrum changes with time.
import numpy as np
import os
import matplotlib as mpl
mpl.use('Agg')
import matplotlib.pyplot as plt
# Parameters
FWHM = 0.4 # Full width half maximum
Emin=-16. # Starting energy
Emax= 4. # -----
DE = 0.01 # Energy intervals
dE=np.arange(Emin,Emax,DE)
E_shift = 0
filename = 'liquid_water-k1-1.pdos'
def plot_pdos(txtfile):
if not os.path.isfile(txtfile):
print("Warning: File does not exist!")
return np.zeros(np.shape(dE)[0])
PDOS = np.asarray(open(txtfile).read().split())
if 'f' in PDOS:
start = 26
Norb = 4+3
elif 'd' in PDOS:
start = 25
Norb = 3+3
elif 'p' in PDOS:
start = 24
Norb = 2+3
else:
start = 23
Norb = 1+3
L = int(PDOS[-Norb])
N = np.where(PDOS=='Projected')[0]-1
for i in range(N.shape[0]-1):
pdos = PDOS[N[i]:N[i+1]]
pdos = pdos[start:]
pdos = np.asarray(np.split(pdos,L))
E = pdos[:,1].astype(np.float)*27.2114 #a.u. -> eV
pdos = pdos[:,3:].astype(float)
if i == 0:
pdos_tot = np.zeros((pdos.shape))
pdos_tot += pdos / N.shape[0]
s, p = convolve(E,pdos)
if i == 0:
yshift = np.max(np.append(s,p))
plt.plot(dE+E_shift,s-i*yshift, c='blue',linestyle='-',linewidth=2.,alpha=0.7)
plt.plot(dE+E_shift,p-i*yshift, c='red',linestyle='-',linewidth=2.,alpha=0.7)
s, p = convolve(E,pdos_tot)
plt.text(-14,-(i+1.5)*yshift,'Sum')
plt.plot(dE+E_shift,s-(i+2)*yshift, c='blue',linestyle='-',linewidth=2.,alpha=0.7)
plt.plot(dE+E_shift,p-(i+2)*yshift, c='red',linestyle='-',linewidth=2.,alpha=0.7)
plt.yticks([])
plt.ylabel('PDOS (arbitrary units)')
plt.xlabel('Binding energy (eV)')
plt.xticks(np.arange(-15,5,5),-np.arange(-15,5,5))
plt.xlim([-16,3])
plt.title('FWHM = ' + str(FWHM) + ' eV')
plt.plot([100,101],[0,0], c='blue',linestyle='-',linewidth=2., label='O s')
plt.plot([100,101],[0,0], c='red',linestyle='-',linewidth=2., label = 'O p')
plt.legend(loc=1)
plt.savefig('pdos.png')
plt.show()
def convolve(E,PDOS):
sigma = 1 / (2 * np.sqrt(2 * np.log(2) ) ) * FWHM
conv = np.zeros((len(dE),2))
for j in range(conv.shape[1]):
W = PDOS[:,j]
for i,e in enumerate(E):
conv[:,j] += np.exp(-(dE-e)**2 / (2 * sigma**2))*W[i]
return conv[:,0], conv[:,1]
plot_pdos(filename)