exercises:2019_conexs_newcastle:ex4
Differences
This shows you the differences between two versions of the page.
Both sides previous revisionPrevious revisionNext revision | Previous revision | ||
exercises:2019_conexs_newcastle:ex4 [2019/09/11 03:38] – [Part 2: MD of liquid water] abussy | exercises:2019_conexs_newcastle:ex4 [2020/08/21 10:15] (current) – external edit 127.0.0.1 | ||
---|---|---|---|
Line 298: | Line 298: | ||
&END TOPOLOGY | &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. 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 " | + | 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. |
+ | < | ||
+ | module load Anaconda3/ | ||
+ | </ | ||
+ | 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 " | ||
<code - pdos_conv.py> | <code - pdos_conv.py> | ||
import numpy as np | import numpy as np | ||
- | import matplotlib.pyplot as plt | ||
import os | import os | ||
+ | import matplotlib as mpl | ||
+ | mpl.use(' | ||
+ | import matplotlib.pyplot as plt | ||
# Parameters | # Parameters | ||
Line 380: | Line 390: | ||
< | < | ||
&MOTION | &MOTION | ||
- | &MD | + | |
ENSEMBLE NVT | ENSEMBLE NVT | ||
STEPS 15 | STEPS 15 | ||
Line 395: | Line 405: | ||
&END MOTION | &END MOTION | ||
</ | </ | ||
- | The parameters in the MD section has been selected to finish in a few minut, but feel free to experiment with different values. To not print too much, also modify the &PDOS section to | + | The parameters in the MD section has been selected to finish in a few minutes, but feel free to experiment with different values |
< | < | ||
&PDOS | &PDOS | ||
Line 404: | Line 414: | ||
&END | &END | ||
&END PDOS | &END PDOS | ||
- | which means that we only print the PDOS every third step, as well as appending the specta | + | </ |
+ | which means that we only print the PDOS every third step, as well as appending the spectra | ||
- | When you have finished the calculation, | + | When you have finished the calculation, |
+ | * liquid_water-1.ener - containing information about energies, temperature, | ||
+ | * 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(' | ||
+ | import matplotlib.pyplot as plt | ||
+ | |||
+ | # Parameters | ||
+ | FWHM = 0.4 # Full width half maximum | ||
+ | Emin=-16. | ||
+ | Emax= 4. # ----- | ||
+ | DE = 0.01 # Energy intervals | ||
+ | dE=np.arange(Emin, | ||
+ | E_shift = 0 | ||
+ | filename = ' | ||
+ | |||
+ | def plot_pdos(txtfile): | ||
+ | |||
+ | if not os.path.isfile(txtfile): | ||
+ | print(" | ||
+ | return np.zeros(np.shape(dE)[0]) | ||
+ | |||
+ | PDOS = np.asarray(open(txtfile).read().split()) | ||
+ | |||
+ | if ' | ||
+ | start = 26 | ||
+ | Norb = 4+3 | ||
+ | elif ' | ||
+ | start = 25 | ||
+ | Norb = 3+3 | ||
+ | elif ' | ||
+ | start = 24 | ||
+ | Norb = 2+3 | ||
+ | else: | ||
+ | start = 23 | ||
+ | Norb = 1+3 | ||
+ | |||
+ | L = int(PDOS[-Norb]) | ||
+ | N = np.where(PDOS==' | ||
+ | |||
+ | for i in range(N.shape[0]-1): | ||
+ | pdos = PDOS[N[i]: | ||
+ | pdos = pdos[start: | ||
+ | pdos = np.asarray(np.split(pdos, | ||
+ | E = pdos[:, | ||
+ | pdos = pdos[:, | ||
+ | |||
+ | if i == 0: | ||
+ | pdos_tot = np.zeros((pdos.shape)) | ||
+ | pdos_tot += pdos / N.shape[0] | ||
+ | |||
+ | s, p = convolve(E, | ||
+ | if i == 0: | ||
+ | yshift = np.max(np.append(s, | ||
+ | plt.plot(dE+E_shift, | ||
+ | plt.plot(dE+E_shift, | ||
+ | |||
+ | s, p = convolve(E, | ||
+ | plt.text(-14, | ||
+ | plt.plot(dE+E_shift, | ||
+ | plt.plot(dE+E_shift, | ||
+ | |||
+ | plt.yticks([]) | ||
+ | plt.ylabel(' | ||
+ | plt.xlabel(' | ||
+ | plt.xticks(np.arange(-15, | ||
+ | plt.xlim([-16, | ||
+ | plt.title(' | ||
+ | |||
+ | plt.plot([100, | ||
+ | plt.plot([100, | ||
+ | plt.legend(loc=1) | ||
+ | |||
+ | plt.savefig(' | ||
+ | plt.show() | ||
+ | |||
+ | def convolve(E, | ||
+ | |||
+ | sigma = 1 / (2 * np.sqrt(2 * np.log(2) ) ) * FWHM | ||
+ | |||
+ | conv = np.zeros((len(dE), | ||
+ | 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) | ||
+ | </ |
exercises/2019_conexs_newcastle/ex4.1568173089.txt.gz · Last modified: 2020/08/21 10:15 (external edit)