Liquid water: PDOS and Born-Oppenheimer MD

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.

Part 1: PDOS and spectra broadening of liquid water

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.

liquid_water.inp
&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:

liquid_water.xyz
   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”.

pdos_conv.py
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)

Part 2: MD of liquid water

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.

pdos_conv_md.py
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)