Practical Session 4: Molecular Dynamics - Max Planck Society

1y ago
1.34 MB
7 Pages
Last View : 1y ago
Last Download : n/a
Upload by : Xander Jaffe

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 file is created but terminal is free for other use.– You can look and analyse the output file 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 accuracy rho 1E-5sc accuracy eev 1E-4sc accuracy etot 1E-6All the useful files 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” script to analyse your run, typing in the terminal:perl aims MD "FHI-aims-output-file" "script-output-file"Do this for both outputs. Now plot the total energy (fifth column of the script output) vs.the time of simulation (first column of the script output) in xmgrace. For a short guide onhow to plot files 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 first 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 file 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 flags).A geometry for NH3 is already provided in the folder exercise 1 (all the files for exercisex are located in the folder named exercise x ). You will notice that this file contains alsospecific 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 (flags 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 effects of the time step size, also in a NVE simulation. Keep the same file of the last exercise, and use the same file (withthe tight convergence settings).The only modification 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 anotherfile. Plot the total energy vs. time of simulation (again using the “aims MD”, first and fifthcolumn) 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 fluctuations 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 "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 flag “MD run NVE 4th order” in your, 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) fluctuate2around 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 specified 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 fictitious 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 fluctuations 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 fluctuating 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 and (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 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” 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 file for this simulation canbe found in the folderexercise 4/harmonic vibs.You should use the same 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 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 The syntax is:56

python get 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 file that contains only positions. You should use thesame 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 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 file in the folder:exercise 4/anharmonic vibs/NVE runCopy your 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 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, you need to preparean input file called ‘’. 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, 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 ”” output from the harmonic analysis into the velocity-autocorrfolder and run the script project 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. the name of the MD outputOtherwise, you will need a file with the same flags describedabove in the folder you are running. A typical call will look like:project 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 file containing a thermalized configuration(at 300K) for the cation H5 O 2. In the folder you find a 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 .

Related Documents:

Business Ready Enhancement Plan for Microsoft Dynamics Customer FAQ Updated January 2011 The Business Ready Enhancement Plan for Microsoft Dynamics is a maintenance plan available to customers of Microsoft Dynamics AX, Microsoft C5, Microsoft Dynamics CRM, Microsoft Dynamics GP, Microsoft Dynamics NAV, Microsoft Dynamics SL, Microsoft Dynamics POS, and Microsoft Dynamics RMS, and

Microsoft Dynamics 365 for Operations on-premises, Microsoft Dynamics NAV, Microsoft Dynamics GP, Microsoft Dynamics SL, Microsoft Dynamics AX 2012 or prior versions, or Microsoft Dynamics CRM 2016 or prior versions. This guide is not intended to influence the choice of Microsoft Dynamics products and services or provide technical specification.

This guide is designed to improve your understanding of how to license Microsoft Dynamics 365, Business edition. This document does not apply to Dynamics 365, Enterprise edition, Microsoft Dynamics NAV, Microsoft Dynamics GP, Microsoft Dynamics SL, Microsoft Dynamics AX 2012, or Microsoft Dynamics CRM 2016 or any other prior version.

Microsoft Dynamics Guide Dynamics GP vs. Dynamics 365 (NAV) vs. Dynamics SL . Dynamics 365 BC -1 Premium User 100/month/user (Subscription) 2000 (On-Premise) . Solomon's application became Dynamics SL. Dynamics SL is geared first and foremost for project-based businesses. This makes SL the

Microsoft Dynamics 365 for Operations on-premises, Microsoft Dynamics NAV, Microsoft Dynamics GP, Microsoft Dynamics SL, Microsoft Dynamics AX 2012 or prior versions, or Microsoft Dynamics CRM 2016 or prior versions. This guide is not intended to influence the choice of Microsoft Dynamics products and services or provide technical specification.

The journal Molecular Biology covers a wide range of problems related to molecular, cell, and computational biology, including genomics, proteomics, bioinformatics, molecular virology and immunology, molecular development biology, and molecular evolution. Molecular Biology publishes reviews, mini-reviews, and experimental and theoretical works .

Jan 31, 2011 · the molecular geometries for each chemical species using VSEPR. Below the picture of each molecule write the name of the geometry (e. g. linear, trigonal planar, etc.). Although you do not need to name the molecular shape for molecules and ions with more than one "central atom", you should be able to indicate the molecular geometryFile Size: 890KBPage Count: 7Explore furtherLab # 13: Molecular Models Quiz- Answer Key - Mr Palermowww.mrpalermo.comAnswer key - CHEMISTRYsiprogram.weebly.comVirtual Molecular Model Kit - Vmols - CheMagicchemagic.orgMolecular Modeling 1 Chem Labchemlab.truman.eduHow to Use a Molecular Model for Learning . - Chemistry Hallchemistryhall.comRecommended to you b

Xiangrun's Molecular sieve Tel:86-533-3037068 Website: Molecular sieve Types 3A Molecular sieve 4A Molecular sieve 5A Molecular sieve 13X Molecular sieve PSA Molecular Sieve Activated zeolite powder 3A Activated zeolite powder 4A Activated zeolite powder 5A