Practical Session 4: Molecular Dynamics - Max Planck Society

Use the same control.in and geometry.in (including the velocities that are in there) thatyou used in exercise 1. You will only have to change the MD related flags, explained below.Applications Run the three different thermostats, for 0.5 ps in total, using a 0.001ps time step. Use thescheduler for each thermostat, to set the temperature to 300 and 800 K. The key words youshould use are shown below, and the MD run flag should be erased from your input. Thegeneral syntax of the thermostat flags is:Exercise 4: Harmonic vs. anharmonic vibrations: the NH3 molecule MD command time NVT thermostat name temperature thermostat parameter You should provide an educated guess for the value of each thermostat’s parameter, followingthe explanation given above and the hints below. If you have time, try to play with thesethermostat parameters afterwards. Nose-hooverMD scheduleMD segment 0.15 NVT nose-hoover 300 QMD segment 0.35 NVT nose-hoover 800 QHint: for setting the mass Q of the thermostat degree of freedom, consider it as an oscillator whose proper frequency should be in the range of the vibronic frequency (phonons)of the atomic system. Considering NH3 has frequencies (see next exercise) between 1000and 3500 cm 1 , derive an expression that links Q with its frequency f and estimate Q. (Acontrol file control.nose-hoover programmed.in with a reasonable value of Q is providedin exercise 3/solution/scheduler)Vibrational spectroscopy is, nowadays, a very important tool for the characterization of molecules.Usually, the frequencies measured experimentally are compared to theoretical calculations in orderto determine the geometry and electronic structure of the molecule. However, the most commonapproach for the ab initio calculations consists of relaxing the geometry of a molecule at 0Kand then performing a vibrational analysis in the the harmonic approximation. There are a fewproblems in this approach, for example the inability of probing all representative conformationsof the measured molecule, which usually finds itself at room temperature, and the inability toinclude effects due to anharmonicities that appear in particularly floppy vibrational modes of themolecules.It is possible to go beyond the harmonic approximation by “brute-force”, but calculating the shapeof the potential energy surface even for very few degrees of freedom is an amazingly demandingtask.One can overcome some of these drawbacks by performing a Molecular Dynamics simulation ofthe system in question. In the framework of Linear Response Theory, one can rewrite the FermiGolden rule by means of the Fourier transform of the dipole moment time correlation function [9]: dt eiωt M (t) · M (0) tI(ω) F (ω) (8)In this formula, I(ω) is the intensity of the vibrations and F (ω) is a quantum corrector factor whichmust be multiplied with the classical line shape in order to reproduce measured amplitudes [10, 11].The angular brackets mean a statistical time average of the auto-correlation of the dipole momentof the molecule. Formula 8 will give all frequencies that are active in the IR range. Therefore,within one MD run, the whole IR spectrum of the molecule can be calculated, since one can choosevarious t 0 to average the dipole auto-correlation over. AndersenMD scheduleMD segment 0.15 NVT andersen 300 νMD segment 0.35 NVT andersen 800 νHint: ν (the unit is ps 1 ) should have a value such that every few ( 10) time steps there isa fair chance (20-40%) the system interacts with the bath. BerendsenA similar relation can be found for the time average of the velocity auto correlation function:VDOS(ω) MD scheduleMD segment 0.15 NVT berendsen 300 τMD segment 0.35 NVT berendsen 800 τHint: Use τ (the unit is ps) equal to few time steps. i 1,N dt eiωt vi (t) · vi (0) t(9)In this case, N is the number of atoms in the molecule and the quantities that are computed areall possible frequencies of vibration of the molecule, not only the ones active in the IR range (dueto selection rules). The advantage is that with this function, one can assign the frequencies toindividual atom’s displacements[9]. While waiting for the simulations to complete, you are challenged to demonstrate Eq. 1. Analyse each thermostat’s output by running again the “aims MD eval.pl” script. Plottemperature (column 2) vs. time steps (column 1) and see how the different thermostatsact. Do you see them trying to change the temperature after 0.15 ps? Is the moleculealready equilibrated at the new temperature when the simulation ends? Plot total energy (column 5) vs time step (column 1). Did you expect such behaviour? Now,for the Nose-Hoover case, plot “Nose-Hoover Hamiltonian” (last column) vs time steps (firstcolumn). What do you observe?Time estimation: 30 min1. Calculating the harmonic vibrationsFor this exercise we will use (again) the NH3 molecule. First we will compute the harmonic vibrational frequencies of the molecule in order tocompare with the anharmonic frequencies. The geometry.in file for this simulation canbe found in the folderexercise 4/harmonic vibs.You should use the same control.in file you were using so far with the tight convergencesettings, leaving away all MD flags and adding the line:relax geometry bfgs 5e-3 Run the aims.vibrations.workshop.mpi.pl script, like you did in Day 1 of these tutorial sessions, and obtain the harmonic vibrational frequencies of NH3 . The frequenciesare outputted in the file basic.vib.out with some other information. In order to filteronly the frequencies, use the script get frequencies.py. The syntax is:56

python get frequencies.py basic.vib.out frequencies.datIn order to plot the frequencies like bars in xmgrace, see the Appendix of this handout.2. Calculating the anharmonic vibrationsNow, we will start the MD simulation. It will be done in two parts.First part: thermalization First thermalize the molecule in the temperature that will be simulated. We will simulate three temperatures, namely: 300K, 450K, and 600K. Each group will only doone temperature, due to the time this simulation takes. Ask the tutors to know whichtemperature you get. These are the steps you should follow:Go to the folder:exercise 4/anharmonic vibs/thermalizationThere you’ll find a geometry.in file that contains only positions. You should use thesame control.in file as in exercise 1 with tight convergence settings and add the line:distributed spline storage .true.This line enhances the performance of the parallelization by enabling the distributionof the splined Hartree potential components across the memory of different MPI-tasks.More details about this flag can be found in the FHI-aims manual, page 112. Theoutcome is that the simulation will be a bit faster.We will use the Andersen thermostat and run 0.4ps of the MD NVT simulation in orderto thermalize the molecule. The MD block of your control.in file should look like this:MD run 0.4 NVT andersen your temperature 200MD time step 0.001MD MB init temperatureThe last flag initializes the velocities of the atoms from a random Maxwell-Boltzmanndistribution at a given temperature. You should put the value of the temperatureassigned to you, e.g., 300.Second part: the NVE run After the thermalization run finishes, look at the output, copy the last geometry withthe corresponding atom’s velocities and create a geometry.in file in the folder:exercise 4/anharmonic vibs/NVE runCopy your control.in file that you used in the thermalization in this folder as well. Wewill now run a 3ps NVE simulation with the initial velocities given by the thermalizationthat was just performed. Now you include a important flag for this calculation: outputdipole. Your MD block should read:MD run 3.0 NVEMD time step 0.001This run takes 50 minutes. Dipole-dipole correlation function: IR spectrum. For the analysis of this simulation, we will use the auto-correlate.py script in order to see the evolution of theauto correlation functions and the spectrum with the time of the run. Besides doingthis analysis at the end of the simulation, you can do it after approximately 1ps, 2psof simulation (or more often if you want). In order to interactively check how far thesimulation has gone, you may write:grep ‘‘Updated atomic structure’’ wc -lwhich gives the number of completed time steps, i.e. the time elapsed in fs (in our casewere the time step is 1 fs). In order to run auto-correlate.py, you need to preparean input file called ‘control.autocorr.in’. The flags are the following:– path string‘string’ is the name of FHI-aims output7– choice string‘string’ is the type of autocorrelation you want to calculate. Accepts either ’velocity’or ’dipole’– sampling interval integer‘integer’ is an integer number that defines the sampling interval (each t 0) tocalculate the autocorrelation. Recommended value is 1.– cutoff ratio float‘float’ is a float number in the interval [0, 1] that defines the ratio of the tail ofthe autocorrelation function you wish to leave out in order to make the fouriertransform. Recommended value is 0.1.– broadening float‘float’ is a float number in units of wavenumbers (cm 1 ) that defines the broadeningof the Gaussian to be convoluted with the Fourier transform. Recommended valueis 3.An example of input file for this script is provided in the folder exercise 4/anharmonic vibs/.You run the script by typing:python auto-correlate.pyThree files will be generated, namely:– autocorr.dat containing 3 columns, the first being time in ps, the second being theautocorrelation function, and the third being the autocorrelation function times awindowing function that makes it go to zero on the edges. The windowing functionis essential for reducing the noise of the Fourier transform.– raw fourier transform.dat containing 2 columns, the first being wavenumbers incm 1 and the second the intensities in arbitrary units– convoluted fourier transform.dat containing 2 columns, the first being wavenumbers in cm 1 and the second the intensities in arbitrary units convoluted with agaussian curve. Run at first the python script with the ‘dipole’ option in the dipole-autocorr folder.Remember to rename the outputs of the script so that they don’t get overwritten.Visualize them in xmgrace and see how the autocorrelation function and the peaks ofthe Fourier transform evolve with time. What does the “long” time (average) correlationindicate? Did you expect the curve not to go to zero? Watch the movie. After extracting the xyz movie of the MD run, via the scriptcreate xyz movie.pl, you may also want to visualize the trajectory, e.g., with vmd.Checking the visualized trajectory, together with the plot of the total-should-be-conservedenergy, is generally an unmatched test one does in order to see whether something (andwhat, if the case) went wrong or unexpectedly. Now that you have the anharmonic and the harmonic vibrational frequencies, try toplot them on top of one another to see the differences. The outputs are in arbitraryunits, therefore you should scale one of the two spectra in order to compare. Highertemperatures should show more anharmonic effects, while low temperatures should becloser to the harmonic result. Which peak shows more anharmonicity? Velocity auto correlation function and projection onto normal modes. After the whole FHI-aims run is finished, run also the script with the ‘velocity’ optionin the velocity-autocorr folder. Copy the ”basic.xyz” output from the harmonic analysis into the velocity-autocorrfolder and run the script project vacf.pl in order to project the velocity autocorrelations onto the harmonic normal modes. The parameters of this script are:1) the xyz containing the normal modes (i.e. basic.xyz)2) the name of the MD outputOtherwise, you will need a control.autocorr.in file with the same flags describedabove in the folder you are running. A typical call will look like:project vacf.pl basic.xyz fhi-aims output 8

You will get as output one intensity file for each stable frequency (i.e. without the 6rigid body modes at 0 frequency) of the molecule. The outputs are named:convoluted fourier transform. fhi-aims out .mode. id mode mode freq. .vacfwhere id mode is just a counter that orders the modes.Plot them all together, e.g. usingxmgrace convoluted fourier transform. fhi-aims output .mode.*Are the normal modes still vibrating with a localized frequency? Which mode(s) is/arepresent in the spectrum of the velocity auto correlation function and not in the spectrumof the dipole auto correlation function?Question for music lovers: can you see the Tartini (difference) third “sound”?Timing: 55min for the simulations , 1h15min total Exercise 5: Harmonic vs. anharmonic vibrations: H5 O 2 and D5 O2 (Extra)If you have extra time and curiosity here is an interesting small project (you may run these overnight, while you have fun in Berlin):In folder exercise extra/ you will find a geometry.in file containing a thermalized configuration(at 300K) for the cation H5 O 2. In the folder you find a control.in with optimized settings (notice you will be running GGAPBE now), change it in order to simulate a NVE MD run of 10ps. In case you wish to tryanother temperature, first a 0.5ps run with a thermostat has to be performed. Do the same for the system D5 O 2 (deuterium instead of hydrogen, mass of deuterium:2.01410 a.u.) Calculate the spectrum of the dipole autocorrelation for both H5 O 2 and D5 O2 and comparewith the harmonic frequencies you find in the folder. Which is the most anharmonic feature? Can you assign the modes, by means of the spectrum of the velocity autocorrelation and itsdecomposition? In the folder exercise extra/solution you find results obtained with the “aims-standard”basis set, with a 20ps long NVE MD simulation, at different temperatures, for both H5 O 2and D5 O 2 . Which differencies can you notice with the data you obtained? In the folder exercise extra/ you find also a paper by Sauer and Döbler, published onChemPhysChem 6, 1706 (2005), in which the same system is analyzed in a similar fashionby means of higher level methods, namely MP2 for the anharmonic spectra via MD. How doDFT-PBE results compare? Note: the published results contain also the quantum correctorfactor of Eq. 8, which takes the form: F (ω) κω 2 . The actual form of κ, which justuniformly scales the spectrum, is uninfluential here.Figure 2: The equilibrated H5 O 2 molecule.910

APlotting files in xmgraceIn order to plot files containing multiple Y columns, follow these steps:1. Open xmgrace and click on Data Import ASCII (Figure 3).2. Find the file that you want to plot and choose ”Load as Block Data” in the dialogue boxthat will open (Figure 4).3. Set X from column 1 and Y from whatever other column you want to plot in the new dialoguebox that will open (Figure 5).4. Press Apply and then press Close on all dialogues.In order to plot bars, useful when you want to plot for example the harmonic frequencies as sticks,in step 3 above, instead of leaving “Set type XY”, change for “Set type BAR” (Figure 6).Properties of the plots can be changed by double clicking on the data or on the axis.Figure 4: Step 2 - Find the file that you want to plot and choose ”Load as Block Data” in thedialogue box that will open.Figure 3: Step 1 - Open xmgrace and click on Data Import ASCIIFigure 5: Step 3 - Set X from column 1 and Y from whatever other column you want to plot inthe new dialogue box that will open.1112

References[1] D. Frenkel and B. Smit, Understanding Molecular Simulation: from algorithms to applications,second edition, Academic Press 2002.[2] A. M. N. Niklasson, C. J. Tymczak, and M. Challacombe, Phys. Rev. Lett. 97, 123001 (2006)[3] T. A. Arias, M. C. Payne, and J. D. Joannopoulos, Phys. Rev. B 45, 1538 (1992)[4] A. M. N. Niklasson, Phys. Rev. Lett. 100, 123004 (2008)[5] H. Ishida et al., Chem. Phys. Lett. 282, 115 (1998)[6] H. C. Andersen, J. Chem. Phys. 72, 2384 (1980)[7] S. Nosé, J. Chem. Phys. 81, 511 (1984). W.G. Hoover, Phys Rev. A 31, 1695 (1985)[8] H. J. C. Berendsen et al., J. Chem. Phys. 81, 3684 (1980)[9] M-P. Gaigeot, M. Martinez, and R. Vuilleumier, Mol. Phys. 105, 2857 (2007)[10] J. Borysow, M. Moraldi, and L. Frommhold, Molec. Phys. 56, 913 (1985)[11] R. Ramirez, T. Lopez-Ciudad, P. Kumar, and D. Marx, J. Chem. Phys. 121, 3973 (2004)Figure 6: Useful to plot harmonic frequencies as sticks.1314

2. Extended Lagrangian: Nose-Hoover thermostat [1, 7]. Equations of motion derived from the Lagrangian of the system conserve the total energy of the system. One can write an extended Lagrangian, by adding fictitious degrees of freedom, such that the overall total energy is conserved but the atomic subsystem can span ensembles .