1y ago

15 Views

1 Downloads

1.34 MB

7 Pages

Transcription

sc accuracy rho 1E-2sc accuracy eev 1E-1Practical session 4: Molecular dynamicsFriday, June 26 2009, 14:00sc accuracy etot 1E-3sc accuracy forces 5E-2In order to start the run type:mpirun "FHI-aims binary" -np 2 "output file" &– The & puts the run on background.– Output ﬁle is created but terminal is free for other use.– You can look and analyse the output ﬁle while the program is running.– ATTENTION: do not start another FHI-aims run simultaneously.It will be VERY slow. When the previous calculation is over, run another simulation, changing the name of theoutput and using:Figure 1: The NH3 molecule.sc accuracy rho 1E-5sc accuracy eev 1E-4sc accuracy etot 1E-6All the useful ﬁles for this tutorial (together with solutions) are located in:/usr/local/aimsfiles/tutorial4We strongly suggest to copy the whole tree into your home directoryand use the copy as working directories:cp -r /usr/local/aimsfiles/tutorial4.All the scripts for this tutorial are located, as usual, in:/usr/local/aimsfiles/scriptssc accuracy forces 5E-4When it is done, use the “aims MD eval.pl” script to analyse your run, typing in the terminal:perl aims MD eval.pl "FHI-aims-output-file" "script-output-file"Do this for both outputs. Now plot the total energy (ﬁfth column of the script output) vs.the time of simulation (ﬁrst column of the script output) in xmgrace. For a short guide onhow to plot ﬁles with multiple columns in xmgrace see the Appendix of this tutorial.Can you see how the energy drifts with the not-so-accurate settings?The NVE ensembleExercise 1: The importance of the SC convergence criteriaTo start our ﬁrst Molecular dynamics simulation, we will use a very simple molecule, which isAmmonia (NH3 ), and we will investigate the importance of the self-consistency convergence criteriawhen simulating the NVE ensemble.Ideally there should be no energy drift whatsoever since in the NVE ensemble, the energy shouldstay constant. The reason for this drift is that we go away from the true Born-Oppenheimer surfaceif we don’t converge well our quantities and this leads to an unphysical (and undesirable) energydrift.Timing: 5 minutes per simulation, 30 minutes total First, build an input ﬁle for FHI-aims using the LDA (pw-lda) functional, no spin polarization, specifying the use of Molecular Dynamics in the NVE ensemble, and using the “light”standards for the species (please refer to the manual for the exact syntax of these ﬂags).A geometry for NH3 is already provided in the folder exercise 1 (all the ﬁles for exercisex are located in the folder named exercise x ). You will notice that this ﬁle contains alsospeciﬁc velocities for each atom, so that the Molecular Dynamics run will not use a randominitialization. These velocities come from a previous equilibration of the molecule at 300K.You will see how this thermalization is done in exercises 3 and 4. Start a 0.15ps NVE run, using a 0.0005ps (0.5fs) time step (ﬂags MD run and MD time steprespectively), with the following self consistency convergence criteria:12

Exercise 2: The importance of the time step sizeIn the second exercise, we will investigate the eﬀects of the time step size, also in a NVE simulation. Keep the same geometry.in ﬁle of the last exercise, and use the same control.in ﬁle (withthe tight convergence settings).The only modiﬁcation to be made is to increase the time stepof the MD simulation to 0.001ps, and run FHI-aims again. Now increase the time step to 0.003ps and run the program, directing the output to anotherﬁle. Plot the total energy vs. time of simulation (again using the “aims MD eval.pl”, ﬁrst and ﬁfthcolumn) for the result you obtained in exercise 1 for 0.0005ps time-step and tight settings,and the other two you just obtained (0.001ps and 0.003ps time steps) in the same plot,again in xmgrace. How do the energy ﬂuctuations develop? What happens with the 0.003pstime-step run?It is interesting to have a bigger time step because one could calculate longer trajectories in a shortertime. Notice, however, that 0.003ps simulation diverges. In fact, if you look at the evolution ofthe geometry, by runningcreate xyz movie.pl "FHI-aims-output-file" "script-output-file".xyzon your output and opening it in VMD, you will see that the molecule dissociates. This happensbecause the integrator is unable to deal with these “big” time steps. This integrator uses a simpleVerlet algorithm [1], where the error in the trajectory goes with Δt4 .Now we will investigate how higher order integrators perform, in particular, the fourth orderintegrator as proposed by Ishida et al.[5]. Use the ﬂag “MD run NVE 4th order” in your control.in, choose a 0.002ps time step, andrun the program again using all same settings as before. Look at the output and see thenumber of force evaluations needed in this case. The line where this information is outputtedlies in the end of the FHI-aims output. It should be the same as using the usual integratorwith a 0.0005ps time step (as was done in exercise 1).The NVT ensembleExercise 3: Testing thermostatsFrom a statistical mechanics point of view, the temperature is imposed by bringing the systeminto thermal contact with a large heat bath. The velocities (linear momenta) of the particles inthe system are distributed according to the Maxwell-Boltzmann (MB) distribution: P (p) β2πm 3/2 exp βp2 /(2m)(1)This means that the Temperature (kinetic energy) is not constant but can (and should) ﬂuctuate2around the average value. The theoretical standard deviation is 3N, where N is the number ofatoms.Here we apply three schemes which simulate a thermostat in MD. Below is a short summary ofthe thermostats we will use in this tutorial:1. Stochastic thermostat: Andersen thermostat [1, 6].It tries to apply literally the concept of “coupling with a heat bath”. In practise, occasionallya particle is randomly selected and its velocity is drawn from the MB distribution at the targettemperature. The algorithm is speciﬁed by the temperature and by the coupling parameterν, such that the probability that a particle is selected in a time step Δt is νΔt.2. Extended Lagrangian: Nose-Hoover thermostat [1, 7].Equations of motion derived from the Lagrangian of the system conserve the total energy ofthe system. One can write an extended Lagrangian, by adding ﬁctitious degrees of freedom,such that the overall total energy is conserved but the atomic subsystem can span ensemblesother than microcanonical. With the Nose-Hoover Lagrangian, the atomic subsystem samplesthe canonical ensemble. The equations of motion of the Nose-Hoover thermostat are:ṙiṗi Plot the total energy vs. time for the 4th order integrator with 2fs time step and the normalverlet integrator with 0.5fs time step. How do the energy ﬂuctuations compare?η̇π̇Timing: 20 minutes total pi /mi U rNπpi riQπ Q p2gi miβi(2)(3)(4)(5)where g is the number of degrees of freedom of the system. The conjugated momentum π ofthe extra coordinate η acts as a ﬂuctuating drag parameter to the atomic subsystem. Theconserved energy associated to the equations of motion is:E p2 1 π2ηi g U rN 2m2Qβii(6)3. Berendsen temperature coupling (not sampling NVT!) [1, 8]According to the Berendsen algorithm, the temperature of the system is controlled by scalingthe velocities with a factor λ: 1/2 Δt T0λ 1 1(7)τTwhere T0 is the target temperature, T is the actual one, Δt is the time step and τ is aparameter that controls the strength of the coupling. If τ Δt the scheme brings back theactual temperature to the target one exactly at each time step. If τ Δt the action of therescaling is milder and milder with increasing τ .34

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 ﬂags, explained below.Applications Run the three diﬀerent 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 ﬂag should be erased from your input. Thegeneral syntax of the thermostat ﬂags 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 ﬁle 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 ﬁnds itself at room temperature, and the inability toinclude eﬀects due to anharmonicities that appear in particularly ﬂoppy 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 diﬀerent 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 (ﬁrstcolumn). 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 ﬁle for this simulation canbe found in the folderexercise 4/harmonic vibs.You should use the same control.in ﬁle you were using so far with the tight convergencesettings, leaving away all MD ﬂags 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 ﬁle basic.vib.out with some other information. In order to ﬁlteronly 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 ﬁnd a geometry.in ﬁle that contains only positions. You should use thesame control.in ﬁle 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 diﬀerent MPI-tasks.More details about this ﬂag 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 ﬁle should look like this:MD run 0.4 NVT andersen your temperature 200MD time step 0.001MD MB init temperatureThe last ﬂag 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 ﬁnishes, look at the output, copy the last geometry withthe corresponding atom’s velocities and create a geometry.in ﬁle in the folder:exercise 4/anharmonic vibs/NVE runCopy your control.in ﬁle 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 ﬂag 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 ﬁle called ‘control.autocorr.in’. The ﬂags 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 deﬁnes the sampling interval (each t 0) tocalculate the autocorrelation. Recommended value is 1.– cutoff ratio ﬂoat‘ﬂoat’ is a ﬂoat number in the interval [0, 1] that deﬁnes 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 ﬂoat‘ﬂoat’ is a ﬂoat number in units of wavenumbers (cm 1 ) that deﬁnes the broadeningof the Gaussian to be convoluted with the Fourier transform. Recommended valueis 3.An example of input ﬁle for this script is provided in the folder exercise 4/anharmonic vibs/.You run the script by typing:python auto-correlate.pyThree ﬁles will be generated, namely:– autocorr.dat containing 3 columns, the ﬁrst 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 ﬁrst being wavenumbers incm 1 and the second the intensities in arbitrary units– convoluted fourier transform.dat containing 2 columns, the ﬁrst being wavenumbers in cm 1 and the second the intensities in arbitrary units convoluted with agaussian curve. Run at ﬁrst 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 diﬀerences. The outputs are in arbitraryunits, therefore you should scale one of the two spectra in order to compare. Highertemperatures should show more anharmonic eﬀects, 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 ﬁnished, 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 ﬁle with the same ﬂags 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 ﬁle 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 (diﬀerence) 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 ﬁnd a geometry.in ﬁle containing a thermalized conﬁguration(at 300K) for the cation H5 O 2. In the folder you ﬁnd 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, ﬁrst 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 ﬁnd 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 ﬁnd results obtained with the “aims-standard”basis set, with a 20ps long NVE MD simulation, at diﬀerent temperatures, for both H5 O 2and D5 O 2 . Which diﬀerencies can you notice with the data you obtained? In the folder exercise extra/ you ﬁnd 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 uninﬂuential here.Figure 2: The equilibrated H5 O 2 molecule.910

APlotting ﬁles in xmgraceIn order to plot ﬁles containing multiple Y columns, follow these steps:1. Open xmgrace and click on Data Import ASCII (Figure 3).2. Find the ﬁle 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 ﬁle 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 ﬁctitious degrees of freedom, such that the overall total energy is conserved but the atomic subsystem can span ensembles .

Related Documents: