Self-Examination Questions for CH3F1 Molecular Modelling PDF

Self-Examination questions for CH3F1 Molecular Modelling Answers are only given for questions that are not direct review of lecture material. 1. Explain the terms theory, model, and reality and clarify their differences? Give examples 2. What is the purpose of MM and how does it contribute to chemis...

Self-Examination questions for CH3F1 Molecular Modelling Answers are only given for questions that are not direct review of lecture material. 1. Explain the terms theory, model, and reality and clarify their differences? Give examples 2. What is the purpose of MM and how does it contribute to chemistry research? 3. Explain the concepts of Inverse Design Molecular Modeling and Multiscale Modeling and how they differ from usual Molecular Modeling. 4. Explain the difference between high performance computing and distributed computing, such as Folding@Home 5. Give a step-by-step explanation of the Born-Oppenheimer approximation and its practical consequences for molecular modeling and interpretation of chemistry. 6. Can you think of cases where the Born-Oppenheimer approximation would not be valid and might break down? Model Answer: The BOA is not valid in cases where the nonadiabatic coupling becomes significantly large. This is the case, when nuclei move very fast due to external forces or when electronic wave functions vary strongly due to nuclear motion. The best example is photochemistry (light-induced chemistry), where atomic motion is induced by electronic excitation. 7. Write down the full Schrödinger equation with explicit terms for the Hydrogen molecule-ion H2+. Explain each term. Model Answer: This expression is identical to the one mentioned in the lecture, with the exception of the fact that now there is only 1 electron, so there is no sum over electrons in the electronic kinetic energy and the electron-nuclear repulsion and, most importantly, there is no electron-electron repulsion term. Note: Try to write down this equation for practice at least once. 8. Explain the vibrational enhancement effect in dissociative adsorption at hand of an ‘elbow’ plot. Model Answer: Please review lectures 2 and 3 9. Explain the basic properties of ab-initio and molecular mechanics methods and clarify their differences and advantages. Model Answer: Ab-initio methods are approximate, but accurate solutions to the electronic Schrödinger equation. They involve complicated mathematic operations (matrix algebra) and quickly become computationally unfeasible for systems larger than 100 atoms. They provide access to electronic and spectroscopic properties of the system and when using them, we do not need to know much about the chemistry beforehand. Molecular mechanics (force fields, interatomic potentials), on the other hand, are empirically-chosen analytical mathematical functions that model the potential energy surface of a molecular system. They are easy to evaluate and compute. Unfortunately, they are very system-specific and we need to know much about the chemistry beforehand. They also involve predefined parameters. 1 CH3F1, 2018/2019 10. Optimization methods use the mathematical properties of PES functions to find minima and 1st order saddle points on a potential energy surface. What are these mathematical properties of minima and 1st order saddle points on a PES? Recall 1st year math! Model Answer: A PES minimum has the lowest energy and is only surrounded by points with higher energy. At the point of a minimum, the force (negative first derivative with respect to the atomic positions) is zero and the curvature (second derivative with respect to the atomic positions) is positive. A 1st order saddle point corresponds to a transition state between a reactant state and a product state and is a point at which the force is zero and the curvature is positive in all but one direction where it is negative. 11. What is the chemical significance of PES minima and saddle points? 12. Sketch the basic ideas behind hybrid QM/MM models. Name possible example applications. 13. Explain when and how to use periodic boundary conditions in molecular modelling? Model Answer: Periodic boundary conditions (PBCs) should be used for systems in which translational symmetry exists, for example, solids, molecular crystals, but also for condensed matter systems such as liquids and amorphous solids where no clear local order exists. Periodic boundary conditions enable the simulations of extended condensed systems without having to describe a large macroscopic number of atoms. In order to use PBCs, we define a unit cell via lattice vectors L and all properties outside this unit cell become images of properties inside the unit cell. Formally P(x) = P(x+L). 14. Describe in detail the terms topology, parametrization, and functional form in the context of a conventional empirical force field for methanol. 15. Explain the optimal choice of topology and functional form of a force field to simulate protein folding dynamics in liquid water in terms of computational efficiency and chemical validity. Model Answer: In protein folding, we are interested in conformational changes of a molecule, therefore we do not need a force field that can make or break bonds (reactive force field). A conventional force field with a topology that contains intramolecular bond stretch, angle bending, and dihedral torsion terms and intermolecular electrostatic and vander-Waals terms will be best choice in terms of efficiency and accuracy. In terms of functional form typical harmonic terms for intramolecular degrees of freedom will be sufficient, as well as a Lennard-Jones-type Van-der-Waals interaction. We also need an efficient and simple water model to enable a simulation of 1000s of atoms to describe a solvated protein. 16. Explain the optimal choice of topology and functional form of a force field to simulate the solid to liquid transition of a metal in terms of computational efficiency and chemical validity. Model Answer: The force field needs to be able to describe isotropic metallic bonding, therefore Many Body potentials or Embedded Atom Models (EAMs) would be the best choice, as they are better at describing highly coordinated bonding scenarios. They also represent a much more efficient approach than actually using ab-initio methods (approximations to the Schrodinger equation). EAMs can model metal structures very 2 CH3F1, 2018/2019 accurately, but in the parametrization care has to be taken to be able to describe properties in the solid and liquid phase. 17. Explain the step-by-step procedure of force field parametrization for an FF that describes liquid water and explain how to ensure its transferability to describe ice. Model Answer: During parametrization, we start from an initial set of parameters that either derive from experiment or ab-initio simulation to best describe individual potential terms. We first calculate target properties by means of molecular simulation and then compare to experiment and ab-initio reference data. We then vary parameters and repeat the calculation and comparison until the agreement with experiment or ab-initio is at a desired accuracy. 18. On the example of a water dimer, write down all terms of a typical empirical force field. Collect all required parameters in a list and break the list down to the set of parameters that need to be fitted. 19. Show mathematically the bond-breaking ability or lack thereof for a harmonic potential and a Morse potential by calculating the force (negative 1st derivative with respect to distance) acting on the bond at infinite atom separation. Model Answer: A bond can break if there are no attractive forces acting between the two atoms at large separation. The harmonic potential is: 𝑉𝑉ℎ𝑎𝑎𝑎𝑎𝑎𝑎 (𝑟𝑟) = 1 𝑘𝑘(𝑟𝑟 − 𝑟𝑟 0 )2 2 Its force is the negative first derivative is (using the chain rule of differentiation): 𝐹𝐹ℎ𝑎𝑎𝑎𝑎𝑎𝑎 = − 𝑑𝑑 𝑉𝑉𝑠𝑠𝑠𝑠𝑠𝑠 (𝑟𝑟) = −𝑘𝑘(𝑟𝑟 − 𝑟𝑟 0 ) 𝑑𝑑𝑑𝑑 This means that at for larger r the attractive force between the atoms becomes more and more negative forcing the atoms closer together. The Morse potential is 0 2 𝑉𝑉𝑀𝑀𝑀𝑀𝑀𝑀𝑀𝑀𝑀𝑀 (𝑟𝑟) = 𝐷𝐷 1 − 𝑒𝑒 −𝛼𝛼 𝑟𝑟−𝑟𝑟 Its force is (using the chain rule and the rule for exponential derivatives): 𝐹𝐹𝑀𝑀𝑀𝑀𝑀𝑀𝑀𝑀𝑀𝑀 (𝑟𝑟) = − 𝑑𝑑 𝑉𝑉𝑀𝑀𝑀𝑀𝑀𝑀𝑀𝑀𝑀𝑀 0 0 = −2𝐷𝐷𝐷𝐷[1 − 𝑒𝑒 −𝛼𝛼 𝑟𝑟−𝑟𝑟 ]𝑒𝑒 −𝛼𝛼(𝑟𝑟−𝑟𝑟 ) 𝑑𝑑𝑑𝑑 This expression gives zero force only at the minimum position 𝑟𝑟 = 𝑟𝑟 0 and at infinite distance 𝑟𝑟 → ∞ when 𝑒𝑒 −2𝛼𝛼 𝑟𝑟−𝑟𝑟 0 goes to zero. This means a bond between two atoms can dissociate. 20. What are the general properties of empirical force fields in contrast to ab-initio methods? 21. Explain the difference in handling intramolecular force field terms and intermolecular force field terms. What additional complications arise for interatomic force field terms? Model Answer: Non-reactive intramolecular force field terms (excluding the Morse potential) are very short ranged and easy to calculate by summing over the topology (list of 3 CH3F1, 2018/2019 interaction terms). Intermolecular terms, such as van-der-Waals potential and the electrostatic interaction are long range and decay slowly. This means that many pairs of interactions need to be included over long distances. Whereas the van-der-Waals potential can be truncated after a certain distance, the electrostatic potential cannot and special algorithms have to be employed that typically are the most time-consuming part of a forcefield calculation. 22. Explain the following: I. What are the general properties of the embedded atom method? II. Describe in words what the idea and functional form behind this specific potential type is. III. Name a material for which this type of potential is particularly suitable 23. Explain the step-by-step procedure of setting up a computer experiment to study the lowtemperature reconstruction of a Au(111) surface. Make your own choices regarding suitable model size and force field method and explain your reasoning. Model Answer: I The Hypothesis or research question is to study reconstruction at a surface. This means we need to find the thermodynamically most favorable, stable surface structures. II Model design: (a) Structure: We have to model the system in a periodic surface slab that is large enough to enable to describe many different reconstructions, ie. the slab has to be much larger than the primitive unit cell. (b) Method: We need to select a method that correctly describes the structural properties of metals and is computationally efficient enough to enable the study of many different structures. From our discussed methods, the Embedded Atom Method would be the most fitting choice. III Simulation: We can either use global optimization methods to search for the energetically most favorable structures or we use Molecular Dynamics simulation to dynamically sample all possible configurations. IV Analysis: We will collect all structures and their relative energies and compare them to data from diffraction and spectroscopy. 24. Explain which assumptions and approximations are taken when simulating atomic motion with classical molecular dynamics and under which conditions they are valid? Model Answer: In classical molecular dynamics (MD), we assume that atomic motion can be described as classical, ie. the atoms follow Newton’s equations-of-motion. The forces that they experience are derived by the potential energy surface that is either calculated via force-fields or ab-initio methods. We thereby assume (a) that the Born-Oppenheimer approximation holds and (b) that quantum effects due to the wave function nature of atoms (nuclear wave function) can be neglected. This includes tunneling effects and vibrational zero point energies. This is typically the case for atoms heavier than Hydrogen and also Hydrogen under certain conditions, but at low temperatures and high vibrational frequencies, this approximation might not be valid. 25. What makes a good integration algorithm and what defines the choice of time step in an MD simulation? 4 CH3F1, 2018/2019 26. What are the consequences of choosing a time step that is too large or too small? 27. We attempt to simulate the protein folding dynamics of an enzyme with 5000 atoms with the Velocity Verlet method. To capture the full dynamics we need to simulate at least 2 milliseconds. The minimal time step we can use is 1 femtosecond (This is limited by the highest vibrational frequency in the system, which typically are C-H stretch vibration.). Every time step takes about 25 microseconds per step of computing time on a simple Laptop. How many MD steps to we have to run and how much Laptop computing time will it take? How long will it take if we can run the simulation on 8 Laptops in parallel? Model Answer: 2 milliseconds = 2 ∙ 10−3 seconds. 1 femtosecond = 1 ∙ 10−15 seconds Therefore, we need 2 ∙ 1012 steps (2 trillion steps). This corresponds to 5 ∙ 1013 microseconds of actual computing time on a Laptop. 1 microsecond = 1 ∙ 10−6 seconds. That means we need 5 ⋅ 107 seconds or ca. 13900 hours or 578 days. Assuming that this could be perfectly parallelized over 8 laptops, this could mean it takes 72 days. Parallelizing over 128 computers, would yield a runtime of 4 ½ days. Note: 25 microsecond evaluation time is very low. Typically this number is larger. Evaluation time is not easy to reduce. The best guess to make simulations faster is rather to increase the time step by choosing better integration algorithms and other tricks. 28. Explain the difference between hot and cold initialization and mention a method to deposit a certain amount of kinetic energy in particles. 29. How do you sample initial conditions for a molecule, which is vibrationally excited but otherwise doesn‘t move? Model Answer: This means that we need to place energy into vibrational motion, but not otherwise induce rotations or translations (at least at first). The easiest way to do this for a diatomic molecule is to look at its bond stretch potential and study the displacement that is necessary in the potential to increase the potential energy to a certain value. We then increase the equilibrium distance of the two atoms by this displacement value. That means we deposited all energy in form of potential energy (momenta are zero, ie. kinetic energy is zero, “cold start”). If we then run an NVE MD simulation of this isolated molecule, we can collect different configurations along the run and these will, after an equilibration time, form a suitable distribution of momenta and positions, which we can use as initial conditions for another simulation. 30. What are the two different stages of a molecular dynamics run and which one is used to extract averaged properties? How can we tell the difference between the two stages? Model Answer: We distinguish between an equilibration phase and a production phase. We can only extract properties from the production phase to generate unbiased time-averages. We can study fluctuations in total energy and temperature to judge if sufficient equilibration time has passed or we can study the change of an averaged value (averaged over short time windows) as a function of time, but the most rigorous way to test if equilibration is sufficient (and if ergodicity can actually be achieved) is to study the velocity-velocity autocorrelation function. 31. Explain the concepts of convergence, ergodicity and statistical sampling. Give examples for non-ergodic systems. 5 CH3F1, 2018/2019 32. Let‘s draw 1000 energies from a production trajectory from MD with 100,000 steps. In which way should we select these energies to minimize the statistical error of the energy average? Explain your choice. Why not use all energies? 33. How much does statistical error decrease by doubling the length of the MD trajectory? 34. How would you calculate the self-diffusion coefficient of pure water? Model Answer: We would calculate the velocity autocorrelation function and integrate it over sufficient time. The area under this function represents the diffusion coefficient. ∞ ∞ 𝑡𝑡=0 𝑡𝑡=0 1 1 𝐷𝐷 = 𝑑𝑑𝑑𝑑 𝐶𝐶𝑣𝑣𝑣𝑣 (𝑡𝑡) = 𝑑𝑑𝑑𝑑 ⟨𝑣𝑣(𝑡𝑡)𝑣𝑣(0)⟩ 3 3 Note: Another correct answer would be to calculate self-diffusion via the averaged meansquared displacement (MSD) of water molecules (so-called Einstein relation), but this was not discussed in the lecture. 35. What is the significance of the velocity autocorrelation function Cvv at t=0? 36. Explain in detail the difference between static and dynamic observables and name examples for both 37. Decide if the following properties can be calculated as dynamic or static observables. Explain your decision: I. Adsorption energy of a small mobile peptide at Au(111) measured at 300K II. Experimentally measured energy loss of a molecular beam after scattering from a surface III. Protein crystal structure at 180 K inferred from Small Angle x-ray scattering IV. Experimentally measured temperature at which molecules desorb from a surface V. Diffusion rate of atomic Hydrogen at a Fuel Cell Cathode VI. Diffusion rate of atomic Hydrogen at an electrode during electrolysis Model Answer: I. Static. The adsorption energy should be independent of time if it is at equilibrium. II. Dynamic. We have to study the actual dynamical motion of the particles. This can be expressed as a reaction correlation function (side-side correlation function) and therefore is dynamic. III. Static. The structure of a molecule or crystal should be independent of time if it is at equilibrium. IV. Dynamic. The desorption temperature is the lowest temperature at which the probability of desorbing from the surface is 100%. We therefore calculate the desorption probability (a dynamic property) as a function of time to identify this temperature. V. Dynamic. We calculate the velocity autocorrelation function. VI. Dynamic. We calculate the velocity autocorrelation function. This is independent of the presence of external electric fields as is the case in an electrochemical cell. 6 CH3F1, 2018/2019 38. Try to think of a method to simulate the lifetime of the intramolecular stretch vibration of a carbonyl group in a large metal-organic complex. (Excitation means that this vibrational mode contains more energy than the other degrees of freedom). Model Answer: The lifetime is the time it takes to loose half the energy to other intramolecular vibrations. If we assume that the vibration loses energy exponentially, this would mean 𝐸𝐸(𝑡𝑡) = 𝐸𝐸0 𝑒𝑒 −𝜏𝜏𝜏𝜏 Here E is the current energy and 𝐸𝐸0 is the initial energy. 𝜏𝜏 is the decay rate. The lifetime is the time at which 𝐸𝐸 = 𝐸𝐸0 2. To determine this we need to determine the decay rate 𝜏𝜏. This is clearly a non-equilibrium dynamic property, we therefore can not use time-averages, but must use ensemble averaging. This means we need to setup an ensemble of initial conditions from which we start the simulations. This is best done by running a long constant energy (NVE) MD simulation from which we then select initial conditions, where every motion in the molecule has the same energy (aka the molecule is in equilibrium). We then take these initial conditions and deposit excess energy into the carbonyl intramolecular stretch (see answer to 29) and use this as starting point for a MD simulation. We do this for many starting points and for each simulation, we collect the total energy of the vibration (sum of kinetic and potential) as a function of time and try to fit it with in exponential decay. We ensemble average over all our starting points and arrive at an exponential decay rate. This we use to calculate the lifetime of the vibration, which we can directly compare with transient spectroscopy measurements. 7 CH3F1, 2018/2019