CP2K Open Source Molecular Dynamics

User Tools

Site Tools

Trace:
exercises:2019_conexs_newcastle:ex4

Differences

This shows you the differences between two versions of the page.

Link to this comparison view

Both sides previous revisionPrevious revision
Next revision		Previous revision
exercises:2019_conexs_newcastle:ex4 [2019/09/11 03:10] – abussyexercises:2019_conexs_newcastle:ex4 [2020/08/21 10:15] (current) – external edit 127.0.0.1
Line 3: Line 3:
 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 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 of liquid water=====+======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.  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.
Line 92: Line 92:
   &END PRINT   &END PRINT
 &END MOTION &END MOTION
 +</code>
 +To run this input, you also need to download the following coordinate file:
 +<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
 +</file>
 +The new section governing the reading of external coordinate files is
 +<code>
 +&TOPOLOGY
 +  COORD_FILE_FORMAT XYZ
 +  COORD_FILE_NAME liquid_water.xyz
 +  CONNECTIVITY OFF
 +&END TOPOLOGY
 +</code>
 +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 
 +<code>
 +module load Anaconda3/5.0.1
 +</code>
 +and then typing
 +<code>
 +python pdos_conv.py
 +</code>
 +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".
 +<code - 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)
 +</code>
 +
 +======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 
 +<code>
 +RUN_TYPE MD
 +</code>
 +in the &GLOBAL section as well as adding the section
 +<code>
 +&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
 +</code>
 +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
 +<code>
 +&PDOS
 +  FILENAME ./liquid_water
 +  APPEND 
 +  &EACH
 +    MD 3
 +  &END
 +&END PDOS
 +</code>
 +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
 +    * liquid_water-1.ener - containing information about energies, temperature, and timings
 +    * liquid_water-pos-1.xyz - containing your trajectory
 +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.
 +<code - 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) 
 </code> </code>
exercises/2019_conexs_newcastle/ex4.1568171439.txt.gz · Last modified: 2020/08/21 10:15 (external edit)