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:22] – [Part 1: PDOS of liquid water] abussy | exercises: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 |
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, | 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, | ||
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 to get pleasing results. | + | 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 | ||
<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 | ||
- | FWHM = 0.3 # Full width half maximum | + | FWHM = 2.5 # Full width half maximum |
Emin=-16. | Emin=-16. | ||
Emax= 4. # ----- | Emax= 4. # ----- | ||
- | DE = 0.01 # Energy intervals | + | DE = 0.01 # Fineness of the energies to convolve to |
dE=np.arange(Emin, | dE=np.arange(Emin, | ||
E_shift = 0 | E_shift = 0 | ||
filename = ' | filename = ' | ||
- | |||
- | ### | ||
def plot_pdos(txtfile): | def plot_pdos(txtfile): | ||
Line 359: | Line 367: | ||
plt.show() | plt.show() | ||
- | | ||
def convolve(E, | def convolve(E, | ||
| | ||
Line 372: | Line 379: | ||
plot_pdos(filename) | 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 | ||
+ | & | ||
+ | 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 ./ | ||
+ | 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, | ||
+ | * 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.1568172140.txt.gz · Last modified: 2020/08/21 10:15 (external edit)