Lecture 32 : Part 2

Another method which is also used called the predictor-corrector algorithm is a high order algorithms that use the transformation or information of the higher order derivatives of the of the atom coordinates. So, instead of taking only the the the d r d t here also we use this let us say or or d 2 d r d 2 r d t square essentially we use del q del q r del r q and all ok. So, the first we do the prediction step. So, we use the predictor. So, the the position velocity and acceleration at time t plus delta t are predicted by the Taylor expansion using using their current value. So, that means, if you want to get t the position velocity and acceleration at time t plus delta t using their position at at at at time t velocity at time t acceleration at at at time t ok. So, these are the kind of the a prediction ok prediction for position prediction for velocity prediction for the for the acceleration. So, as you can see that we are using also some of the higher order differential also ok. Now, higher order differential gives you are using in a prediction case. Then what we are doing that we are calculating the error ok. So, the force is computed at time t plus delta t taking the gradients of the we know the potential at at any position. So, we can do the differential of that potential with respect to position we get the we get the the force ok at t plus t plus delta t and from the force we can easily calculate the the the actual acceleration ok. So, the that resulting acceleration a t plus delta t from this force will be in general difference from the predicted acceleration what is predicted acceleration is here ok and this difference. So, this this is the this is getting we are get from the from the from the negative differential of force and this is from the prediction ok. So, this take the difference of it. So, we get the the error in the in the accelerations. Now, in the second step in the corrected step we use this delta a or this error in measuring or estimating the the the acceleration ok. In the corrected step the positional velocities are corrected proportionally to the calculated errors. So, corrected positions will be the the the predicted position plus a constant into the error in the in the in the in the in the acceleration ok. Similarly, the corrected velocity will be the predicted velocity plus a constant into the of course, this C 0 C 0 and C i will be will be there in the the units will be different to make it a kind of a position and and velocity ok. This constant C i C i means C 0 and C 1 of course, at different atoms will be having this different the different C 0 and C 1 value at a at any step ok. Depend mainly on how many derivatives are included in the Taylor expansion and they generally vary between 1 to 0. This method is very accurate and stable with almost no fluctuation during the run and is especially good for constant temperature molecular dynamics. However, it is not time reversal. Time reversal means whether if you just change plus delta T to minus delta T which whether it is it is it is it is reversible or not ok or remaining same or not ok. However, it is it is not time reversible due to the involvement of the the error correction ok. The, but the earlier methods the the the the valid methods are time reversal. The methods has a tendency to to to show energy drifts with longer time step may be more than 5 femtoseconds and need much more storage space. And this algorithm is further advanced to some of the gear algorithm for higher accuracy by taking higher order derivatives of of atom position. So, these are the three main algorithm which are which are which are used, but apart from this also some of the new algorithms are also coming up ok. Now, knowing all this mainprinciple behind this molecular dynamics, now get to the actual simulation for the how the actual simulation of molecular dynamic done ok. So, essentially start with some of the initial position, initial velocity and structure. Then we go for the initialization ok and also to before we run it also we need to consider some of the the the the pre setup ok. We need to we need to use that the periodic boundary condition. We also need for a long range force calculation also we cannot use let us say force on a particular atom may not be may be dependent on a large number, but we cannot consider all those atoms. So, so we need to use some of the cut off ok and also we need to use some of the nearest neighbor, number of nearest neighbor distance also we need to need to consider the time step ok. Generally, in the most of the cases it is in the range of femtosecond. At the same time also we need to check for the time the the duration for this the the run and some part of the initial run are called the called the equilibration run and then we got the the exact the the the the main MD MD run main MD run from the main MD run only we calculate the the different property ok. So, these are the initialization and some of the precondition then go for the this the the numerical solution of of the numerical also we need to see what are the different external condition or ensemble whether it is a at constant temperature, constant pressure or whether it is at constant energy and all ok and for that essentially we need to use some of the ensembles ok. So, then comes the integration or equilibration step where we use essentially we go for the force calculation from the force calculation we get the the acceleration and also we use those different algorithms ok then we get the data and know this data from the data we can get the difference static and dynamic. So, this is essentially the flowchart for molecular dynamics simulation ok or now as I have mentioned this kind of the parametric or analytical form of the parameterthe potential parameter, but the potential parameter also can be can be generated from the from a from the ab initio calculation. So, those kind of it is this is these are called the ab initio molecular dynamics in that case you do not have to use any external potential external this inter atomic potential, but we can use the the the potential which is generated by the DFT methods and and those potentials can be used to calculate the force and all the using this different algorithm also we can do this the molecular dynamics. So, molecular dynamics are two types essentially classical molecular dynamics where we externally supplied the potential or the molecular dynamics also can be ab initio potential dynamics in that case we need do not have to supply external potential and the potentialgenerally we get the potential from the from the ab initio run directly ok. So, today we are considering mostly the classical molecular dynamics. So, these are the steps for the for the the classical molecular dynamics. So, first we calculate the force on all the atoms from the given potential calculate the the acceleration a on all atoms from the from the from from the force ok and then using this from the. So, if you know the the force from there you get the acceleration and you can use the acceleration to get the velocity and also to get the the new position. So, you can calculate the r i v i and a i. So, a i may be let us say they a kind of it they denote say the position of a particular atom i ok and later time t plus delta t numerically using this finite difference methods and using the data as input for the next round repeat the process under until equilibrium is equilibrium means when the the atoms do not do not subject to a to a to much force ok essentially we consider a kind of a cut off to the force if we say that the force is less than that cut off force we consider that this is under equilibrium ok . So, when the equilibrium is reached and stable trajectories of the atom we get it and from there we calculate the otherproperties ok. The three three most popular numerical integration algorithm we have already discussed or the valid algorithm velocity valid algorithm and also particular algorithm, but nowadays different new methods are also come like gear method and all ok. So, some of the pre setups before we do for this because we need to how many number of atom we need to take in order to represent the true or true materialsystem ok. So, for that we need to go for the some of the the MD pre pre setup one of them is called the periodic boundary condition ok. So, if you consider a simulationyou take a box and consisting of let us say 1000 1000 atoms and let us say they are in a FCC lattice . So, what we can see that we will say that half of this the atom almost 500 atom will be on the surface ok. If it is on the surface then it is coordinated is not same as in the bulk. So, essentially we are not representing the bulk material and all ok. So, how to to to to be considered as a as a bulk atoms of roughly 1 mole you need to. So, in order to consider a bulk you need to have a large number of atom ok. So, that you can neglect some of or you can that means, the number of atoms inside is more than the number of atoms in the in the surface ok. So, in that case you may need may be one Avogadro number of of atoms ok which is kind of impossible toto to deal with ok. So, how to how to how to tackle this problem ok. So, that in order to tackle this problem we use the called the periodic boundary condition take care for this problem by surrounding the box with an infinite number of its image boxes all around. So, you consider a box and then in all the all the faces let us say if you have a cube. So, all these 6 faces even all the on the corner side also it will be all surrounded by some of the imaginary the box we we consider infinite number of imaginary box ok. Then, but the actual simulation takes place only in the atom inside the primary box at the centre centre and all the image boxes just copy the primary boxes ok. For example, for examples so, essentially it is look something like this let us say you are dealing with the let us say this this this central box within using and surrounding the central box we have got a infinite number of similar type of box which are imaginary ok. So, you are running this, but if a particle moves from here to here ok. So, essentially what essentially does that a another particle also move from let us say from the from the below imaginarybelow boxes to the to the centre use in periodic boundary condition. So, actual simulation takes place only for the atom inside the primary box at the centre and all all image boxes just copy the primary box. For example, if an atom k in figures move out of the primary box to the upper image box its image atom in the bottom image box moves in to replace it and keep the number of atom constant in the primary box as I have shown ok. Note that the forces between an atoms and its periodic image can cancel each other and simple coordinate transformation can update the position of atoms in neighboring image boxes if necessary ok. The same trick can be applied not only to the bulk, but also to the single atom and also cluster system by use of this pair. So, not only the the bulk also we canwe can represent asurface or interface by using the periodic boundary condition like here. So, this is the bulk ok. So, all the atoms are essentially we are considering this simple periodicPVC periodic boundary condition ok, but if you have a surface in the case of the surface what essentially the surface is essentially you have got a gap ok. So, what we can do that we can take the periodic boundary conditionperiodicthe the the periodic boundary condition something like this. So, with a periodic box is something like this ok. Similarly, for the case of interface ok let us say we can we can take these two, two different one box one materials another box another materials and all that we can do this ok. So, these are the different representation of PVC periodic boundary condition for or or periodic box for bulk surface and interface ok. So, another piece of also now we now we have got the problem sorted out by using periodic boundary condition to to treat a bulk mostly bulk and surface and all ok. Now, how many. So, in order to calculate the force on a particular atom how many neighboring atom we should consider ok or how many they are the force contribution coming from the neighboring atom also you need to consider. So, for that case essentially we go for the potential cut off ok . The most time consuming part in a molecular dynamics is the is the composition of forces generated by the potential acting on the on the atom essentially the negative differential of potential the gradient negative gradient of the potential . Normally, a potential tail off at long distance as you can see that let us say this is the potential energy curve. So, this at as you as you go go further. So, essentially this you see that the potential energy it is going to going to 0 ok. So, no point of considering a a larger interaction coming from the larger distance ok . So, normally potential tail off at larger distance become negligible and this allow us to disregard the the force calculation where the interactivity distance are beyond a certain cut off distance which call which which we call it the R cut ok. This R cut must be smaller than half of the primary box size it should be. So, this is a convention generally followed it should be half of the primary box size ok . Otherwise an atom may interact with an atom in the primary box and its image in the in the image box which causes some unphysical double counting or double interaction . After the after the cut off the potential tail is normally tapered such that the such that the potential tail vanishes at the at at the at the R cut to prevent any sudden discontinuity on energy and force. Essentially what we are seeing that at R cut the the energy is 0 potential energy is 0 ok , but actually if you go little bit down little bit shorter distance. So, essentially this this is the the energy. So, we make a kind of this tapered case ok . So, this the. So, in this case now the interaction the interactivity force depend only on nearby atoms and overall composition loads scale down to the order of number of atoms ok . So, there are different normal cut offs if you use Lennard-Zurr potential this cut off generally 2.5 to 3.2 into the into sigma . Let take the effective atomic potential or if if atomic methods essentially if 5 angstrom and transfer potential around 3 to 5 angstrom these are the cut offs ok . Also we need to use the minima energy cut off ok . So, essentially let us say . So, this will be cut off let us say this this R cut. So, this this black sphere will be of of of a sphere of radius R cut and this and also first we are considering a periodic boundary condition, then also one atom we are considering whatever the atoms are there within this this this sphere only we consider not the not the other one ok , but at the same time also we need to look for the the the the minima energy criteria also ok. So, during an MD run the force calculation of an atoms are carried out until all the atoms within this potential potential cut off ok. The atom I has a let us say this this atom I has a neighboring atom let us say K ok, but this K within the the main main main calculation box it is outside the R cut off, but if you take the image of this K here. So, you can you can say that this is within the R cut off ok. So, only R cut off is not sufficient because if you when you use the periodic boundary condition what you need to have is called the nearing neighbor neighbor list that also we need to consider because one atom in its main box may be outside the R cut off, but that say the the similar atom in the image on the on the nearby image boxes may be within the R cut off ok. So, that also you need to consider. So, that is what it is mentioned ok . The atom I has a neighboring atom K in the primary box and many K dashed ok like here may be may be somewhere here ok . Above the atom only the atom K dash is in a image box at the at the left in the closest atom closest one to atom I and will be consider for force calculation which is called the minimum minimum image. So, essentially also you need to consider the minimum image criterion . Thus two atoms in the image box at the left are included for force calculation instead of their real real counterparts ok. So, when a potential reaches long range as in the case of the global interaction generated by charge or polarized atoms. So, this criterion is is not not applicable ok . Then also we need to look for the neighbor list. So, during an MV run each atom moves various distances and it is obvious that force calculation for the atoms that moves trivial distances are just meaning meaningless ok . So, Verlet propose a scheme to avoid this unnecessary calculation by listing nearby pairs of atoms and checking their travel distance. And in this technique a larger R circles sphere in 3D of radius R list then the potential cut off is drawn around this Z atom essentially this is the R list ok this one this this gray gray part ok . And in most MD code R list is set around around 10 percent more than the than the R cut as the referred value and a list is made of all the neighbors of each atoms within the within the in the circles. So, this will just avoid this kind of where you can you should not miss this image one ok . After some fixed number of time step for updates out only for the atoms move more than delta R by 2 what what is delta R delta delta R is nothing, but this one this distance this R list minus R cut ok. So, that is the the delta R ok . And it is clear that the list is updated at every every number of this update step before any R listed atoms enter into the into the circle as typically may be in the 10 to 15 molecular dynamic steps ok. So, R list minus R cut essentially more than the number of this this this update number into velocity into delta t ok . Then time step what is the time step you should follow since atoms in a solid lattice vibrates in the range of 10 to the power of minus 14 14 second we have to use minimum time interval in the range of 10 to the power minus 15 because the atomic vibration in the range of 10 to the power minus 14 . So, in order to realistically model the materials the time step has to be less than 10 to the power 14. So, that is why generally we use 10 to the power minus 15 second or 15 to second . And this small time when delta t is assumed that the velocity acceleration and force of each atoms are constant in this time step. So, that we can perform the simple algebraic calculation. So, within this time step there is no change in the acceleration velocity at all ok . To speed up the simulations simulation time step should be increased as much as possible if you increase the time step. So, that to get certain certain times certain time molecular dynamic calculation number of steps should be less ok . So, speed up can be possible by increasing time step, but if you increase time step essential essential things we are also neglect ok . So, a practical rule is that atoms should not travel more than 130 of the nearest neighbor distance in this this time step ok which leads to a typical time step from 0.12 to several femtoseconds ok. Computational load and accuracy with respect to time step is a tradeoff problem ok . So, a smaller time step increase computational load and normally improve accuracy and a larger time step occur the opposite way . If the total energy becomes unstable drift or fluctuates too much it is an indication of too large a time step that overshoot the equilibrium position ok. So, for atoms that move relatively fast due to high temperature or lighter mass or steep potential curve a shorter time step is a common ok. So, time step is also very important generally we use 0.1 to 1 femtosecond ok in molecular dynamics . So, how should be the total time . So, a typical molecular dynamic run may take around 10 to the power 3 to 10 to the power 6 time steps that cause may be million time step you need to go for to total simulation type. So, essentially 10 to the power minus 15 into 10 to the power minus 6. So, essentially we get around 10 to the power minus 9 9 second or some of the nanosecond run. So, that means, we are we are capturing the physics for around 10 to the power minus 9 9 second or nanosecond ok. Normally this is sufficient to see usual static and dynamic properties in materials ok. Running an MD too long may cause error accumulation and inefficiency of data production. Total simulation time however, should be longer than the full relaxation time of the system to have reliable data data data production ok. So, first we initial position or initial velocity are just not not up to the mark or not not at accurate relative to the to the overall atmosphere or overall imposing external condition ok. So, first thing that we need to go for equilibration or relaxation ok. So, that it gives a kind of a relaxed structure and then we run actual run run ok. So, especially for phenomena such as phase transition, vapour deposition, crystal growth and annealing of of irradiation damages the equilibration takes place very slowly and it is necessary to ensure that the total simulation time was long enough. MD run start to bring the unrelaxed initial system to the equilibrium. So, that minimum potential energy where the net force on each atom become 0 under the given condition and so that erasing the memory of the initial. So, first thing that we erase the initial memory initial configuration. So, for that essentially go for first the equilibration step or relaxation step and then the actual the the MD run from there we calculate the properties ok. So, now, in order to take care of the the external condition we need to use the different kind of the ensemble ok. So, in a simulation blocks of a MD run each atom moves and behave differently and a new microstate of the system is generated in every time steps. So, after a proper equilibration and simulation it becomes an ensemble. What is ensemble? It is nothing, but a collection of possible configuration and it has the same macroscopic or thermodynamics property ok. So, in experiment we impose external constraint or ensemble to have the particular property this is especially thermodynamic property such as entropy, Helmholtz free energy, Gibbs free energy and chemical potential can be derived from data obtained from this free energy ensemble. So, in order to get a different particular property we need to use their the their corresponding ensemble ok. So, there are mainly four types of ensembles are there which we call micro canonical called NVE where N is the number of atom, V is the volume, E is the energy is fixed we are not changing it ok. And what kind of system it is? It is a isolated system which is very common in molecular dynamics and in this case we can get the entropy value entropy value from this the partition function if you do this or the ensemble ok. Canonical NVT where energy is not conserved time temperature is conserved energy can vary and this is the most common in in Monte Carlo and common in molecular dynamics we can get the Helmholtz free energy if we run the canonical. Isobaric and isothermal where the the pressure and temperature are are fixed and and also of course, the the number of atoms is fixed and here in this case we can get the Gibbs free energy. Grand canonical in this case the number of atoms may change and and we we can get the then we can get the the chemical potential value ok. So, these are the different ensemble we need to use to get particular type of type of properties ok. And accordingly also algorithm also should be there to to to for a micro canonical system to conserve the the energy for a canonical system it should be some kind of treatment it should be there. So, that the system is at constant temperature because let us say if you temperature is related to to velocity and velocity also vice versa related to to temperature at the same time also velocity are related to pressure also. But pressure we can also change it if we change the change the overall volume also. So, there we need to use different different algorithm to in order to control the temperature. So, this this controlling are called the thermostat or or or or biostat ok. So, these are the different ensemble what we have just described micro canonical canonical isothermal isobaric and and grand canonical. So, except micro canonical essentially what we need to have there you need to have a external bath to conserve the temperature or or to supply of the sufficient particles ok. So, these are the control methods how to control the temperature and and pressure for for for canonical and isobaric and iso isothermal process. So, in an actual experiment a MD run of N-V ensemble often requires a constant temperature and it can be and it can be achieved by rescaling the velocities following the the basic knowledge of kinetic energy ok. So, while kinetic energy we know that is related to the to the velocities or the the temperature ok. And from the average velocity we can we can we are in the temperature also can be directly related to the to the average velocity like this. The average velocity this is kind of N-V ensemble average velocity is over to the square root of temperature in absolute scale ok. So, it is possible to increase the increase or decrease the temperature from T to T dash by multiplying each velocity component. So, from these two you can easily get it ok. Using this equation the velocities of all atoms are gradually increased to the desired value in predetermined iteration steps. And to maintain constant pressure extended ensemble like born states is adopted and and volume of the system is generally allowed to change. So, these are the different the the thermostat and bolster used. Nose hoover, Berenstain, CSVR, Langevin and directly scaling dissipative particle dynamics are popular thermostat and bolster. So, these are mainly nose hoover and Berenstain are used for the for the the for the bolster and others also used for both the mostly for including nose hoover and Berenstain are used for the thermostat ok. So, initialization that means, you need to start from initial structures at the same time also better to give some of the initial velocity also ok. So, student MD to start an MD run several initial input should be prepared to to maintain meaningless to minimize the meaningless data and to have the reliable results ok. So, number of atoms also it should be system size the number of atoms in the primary box should be as small as possible, but at the same time sorry as long as it can represent the real system properly. For reliable data production more than 100 atoms are normally normally required ok. So, initial position and velocity to solve the Newton's equation of motion the initial position and the velocity have to be provided the initial position of the atoms are normally specified according to the known lattice position. The initial velocity of atom are generally normally and randomly chosen from a Maxwell Maxwell Boltzmann or Gaussian distribution at a given temperature. For example, the probability that an atom has a velocity V at a absolute temperature T is given by this expression ok. So, this is the probability that a particular have the velocity V given by m by 2 pi k B T where k B is the Boltzmann constant to the to the power of half exponential minus m V square by 2 k B T ok. The direction of velocities are also chosen randomly to make the the the total linear momentum become 0 ok. So, now, some of the energy minimization steps essentially we use the the numerical method like Steepest Descent method or or Conjugate Variance method ok. In the Steepest Descent method the method start by calculating the potential energy for the initial configuration of atoms and the atoms move towards the lower potential energy until the search line meets the potential gradient perpendicularly. So, you check whether the potential gradient is meeting the next next layer as a perpendicular position ok. Since any further advance leads only to higher higher higher potential energy the next search start at that point along the direction of the Steepest Descent as shown here ok. The process repeats iteratively until the minimum the the configuration minimum energy configuration of the atoms is found. Thus the search path is an zigzag fashion with many orthogonal step that make this this method relatively slow ok. And Conjugate method this method start by searching along the direction of the Steepest Descent and stop at a local minimum and then figure out a new search direction based on both the both the new energy gradients negative value of the force between the atoms and also the previous essentially it is look for the change in the in the gradient. Thus as shown in the the the search path is rather straight forward with fewer step that makes this method relatively faster than the the earlier methods. An extended scheme is the Quasi-Newton method that consider an additional information like second differential of the of the repulsion energy. Now, problem with this energy minimation is that instead of going to a global minima it can also reach to the local minima. That means, this it can have different different localized minima ok and then the local the energy of the localized minima may be more than the global minima. So, you may end up in the local case ok that means, not at the exact very stable position you are reaching to a some kind of the transition state ok. So, in order to avoid that people use different other algorithm it is called the simulated annealing ok. So, after doing all this calculation we get the data and this data can be use of to get the properties. What are the property we can get? We can get the total energy of the of the system given by the summation of the kinetic energy and potential energy at the same time also we get the the the temperature and we can get the the the temperature and energy energy function which and from there from the from the change in the slope we can get the different transition also ok. We can get the structural property like we can get the equilibrium lattice constant, coercive energy, we can get the get the bulk modulus because we are getting the the position, velocity and acceleration and the same time also energy in each step. Maybe we not capturing all the this this this this this this data for a for each step maybe we can we can we can take the average between some of the combined step ok and that we can use for getting this different properties ok. We can also get the thermal expansion coefficient, we can get the radial distribution function. So, this radial distribution function one of the major output from the molecular dynamics are the radial distribution which give the structural characteristics of the of the the the materials ok. So, in the radial distribution as in the X-ray diffraction the radial distribution function here can show how atoms are distributed around a certain atom. So, it is obtained by summing the number of atoms found at a given distance in all direction from a particular atom. So, it just go on increasing the number the the the radius from a particular atom and and count how many number of atoms are there on each incremental and the spherical sphere if we just increase the ok. So, that is give you the radial distribution function and from the radial distribution function we can get what the number of nearest neighbor and also the nearest neighbor distance from the nearest neighbor distance we can say that what kind of bonding it can have and all ok. At the same time also with temperature and pressure how this nearest neighbor the how this this this radial distribution function or r d r d f changes from there you can say that whether any phase change is occurring or not ok like here. Let us say this is the radial distribution functions . So, let us say you have gone increasing the for a particular you take a particular atom ok and then then in c as you as as as you go up away from the there that that that that atom what are the different types the what are the number number of atoms are there ok. So, in this case it is a when you come around let us say 2.2 angstrom. So, you say that number of atoms are count number increases to 5. So, you can say that this is the first coordination ok this is the first this second or you can say the first first nearest neighbor and this is the first neighbor and this is the second neighbor ok . Now, if it is crystalline this will be fixed, but if it is a kind of a amorphous you will say that you get a kind of this sprayed one or if you get a liquid also you will get the this sprayed is more. So, from this radial distribution function we can we can we can say about the the structure of the the material ok the coordination number also if you know the the the the the r d f ok . Also it apart from the r d f also you can get the mean square displacement which can be used to get the the diffusion coefficient and all ok what is mean squared displacement . So, the vector distance that an atom travel may be close to 0 let us let us say it is going left left right up and down, but at the same at a at a time interval it may after if you take the average of it it may become 0. That means, it will say that the average it is it is not moving on a average basis, but at the same time, but molecular dynamics all the atoms are moving ok. So, how to so, you get the for that essentially what we can make it we can make all the displacement all these change in the position as a positive and for that what we take that is a modulus . So, the vector distance that a atom travels may be close to 0 since that atom might take either a forward direction or backward direction or upward direction downward direction ok. Nevertheless if you take the square of the displacement we will have a non zero positive number that is called the mean square displacement or MSD over some time interval t and by arranging the overall time the MSD will have a linear dependence something like this ok we take the square of it . Certain change in the MSD with time are indicative of the melting solidification transition etcetera. In the case of solid the MSD becomes a finite finite value and for a liquid however, atoms will move in indefinitely and the MSD continues to increase linearly with time. The self diffusion coefficient in liquids D can be derived from the slope of the the MSD the time time relationship like this way. What we get the self diffusion coefficient 1 6 of the MSD divided by the the time ok . So, these are the common molecular dynamics simulation programs like GROMACS, LAMBS, NAMD, OpenMM, CHARM, EMBER and all where which mostly use the classical type. That means, where external the the the parametric based analytical the the the empirical polynomials or maybe some functions are used to get the potential to to give the potential parameter ok or molecular dynamics can be of also the where you do not have to supply which is called the ab initio molecular dynamics or AMD some of the program for AMD are CP2K, quantum expresso, VS, VAS, Gaussian, NWK, SIESTA and all ok. So, you may search on this different different program their websites most of them are free ok they are they have also well maintained user manuals and all ok . So, in order to run it essentially we need to go through the their website download it install it at the same time also need to use some of the examples and also go through their the manuals ok . So, how whether the the runs are kind of the good or bad. So, you need to validate it. So, the validation one validation is that conservation of energy ok. So, that is the one way another is that whether confirmation of the global minimum ok and in order to avoid in order to have in order to avoid the local minimum what we can use we can use the simulated annealing. So, the simulated annealing what it is does by equilibrate at higher temperature and slowly cool down to 0 K and the system has enough chances to look look down and identify many different minima on the potential surface and to move over local barriers to arrive at the the global minimum ok . Another validation is the time average under the ergodic hypothesis what essentially say that time and ensemble or statistical average are interchangeable. So, in that case only we can use the statistical method to get the average property ok . So, that is essentially the ergodic you need to confirm to the ergodic or ergodic hypothesis ok and the last one building up or it is diverging and all ok. So, another two types of errors are there in numerical solution one is the transaction error which depends upon the the method itself. So, you need to check you need to change the different methods to to minimize the error or round off error ok. Round off error nowadays is not not much of concern because we are using high precision and the decompression machine and all ok . So, that is conclude that is give the overall view of the principle behind the molecular dynamics simulation there are different flavors of this. Nowadays not only the atoms has a minimum basic unit people are using the number of a a cluster of atom as a unit which is called the coarse grained methods that is also there and there are these molecular dynamics methods not only are used for the materials design or materials property calculation these are also used for in the in the biomolecules medical medical in the industry medicine industry like in the case of the drug design and all ok. Thank you.