 In the previous lecture, we started a little bit of discussion on various aspects of simulation in general. We did not enter specifically into the issues of molecular simulation, but before that we would like to understand that what are the important issues in flow modeling. And one of the very important issues that is of concern is scale issue. So depending on what scale you are addressing, scale means both time scale as well as length scale, different modeling strategies need to be considered. So the fluid modeling or the fluid flow modeling can be broadly classified as continuum modeling and molecular modeling. The continuum modeling for fluid mechanics essentially will usually use the Navier-Stokes equation, but for inviscid flow it will use the Euler's equation of motion. And there are equations which are higher order accurate than the Navier-Stokes equation like for example one order more accurate than the Navier-Stokes equation these are known as Barnett equation. So there are possibilities of using the Navier-Stokes equations or other variants of the same. The molecular modeling on the other hand can be deterministic or statistical. So the deterministic molecular modeling which we will focus on for this particular lecture is molecular dynamics and statistical approaches are also there and two well known statistical approaches are based on like the Boltzmann paradigm and direct simulation Monte Carlo which these are statistical tools. So which modeling strategy should one use? And ideally there is no rigid answer to this particular question, it depends on several other questions. For example what is the system scale, what is the molecular scale, what is the mesoscopic scale that of course everybody knows what is the molecular scale. But question is by setting a question what is the molecular scale, what we essentially intend to mean is that is the molecular scale is capturing the phenomena over the molecular scale important to capture the physics of the problem. So if that is so then one has to invoke molecular modeling or there are some issues which are related to concepts known as multi scale modeling and I will come I will discuss a little bit about those issues if we have some time after we discuss about molecular dynamics. So question is that we generally intend to go for molecular simulations be it molecular dynamics or statistically based molecular simulations if the continuum simulations do not work. Now we have seen and this is just a revisit that we have seen that for gas flows the continuum simulations the continuum hypothesis may not work for high values of the Knudsen number. So typically for Knudsen numbers below 0.1 you can use the continuum modeling safely maybe you may have to adjust the boundary condition with a slip based paradigm instead of a no slip based paradigm but still you can use the continuum considerations with adjustment of boundary conditions. However if you go for higher values of Knudsen number then the continuum nature of the flow is destroyed and discreteness of the molecular entities become important and when this discreteness of the molecular entities becomes important then you have to use typically either molecular dynamics or a Monte Carlo based simulation which is known as direct simulation Monte Carlo or DSMC. So to answer the question whether molecular modeling or continuum modeling will work one has to look into the continuum hypothesis that there should be sufficiently large number of molecules in the chosen elemental volumes. So that if you have uncertainties with regard to their respective positions and velocities that kind of uncertainties do not result in any appreciable deviation in the prediction of flow properties and fluid properties. However the number of molecules present in devices of nanometer dimensions is often lower than that critical value which could have classified the state to be defined by continuum laws. Then the transition from the continuum paradigm to molecular paradigm is decided by a combination of the system length scale and the molecular length scale. For gases the parameter is very easily expressed in terms of the Knudsen number. For liquids molecular mean free path is not a very meaningful parameter. So for liquids the molecular mean free path can be replaced by the parameter sigma which is a finite distance at which inter particle potential becomes 0 and for solid it can be considered to be the lattice constant. With increase of this ratio one needs to employ modifications in continuum theories and finally one has to go for molecular approach. So for gas flows for example we have already discussed about various Knudsen number based regimes and this is just a demonstration that based on the Knudsen number based regime what should be your strategy for simulation. Let us say you are asked to solve a problem of gas flow. Now what kind of simulation will you use? So based on the Knudsen number there is some guideline. So when Knudsen number tends to 0 then you can it is a situation of continuum without molecular diffusion that is without viscous effects. So you can use the Euler's equation less than 10 to the power minus 3 tip these are again as I told these are not hard and fast numbers just typical orders. You can use Navier-Stokes equation with no slip boundary condition at the wall. Between 10 to the power minus 3 to 0.1 you can use Navier-Stokes equation with first order slip boundary conditions at the wall. The first order slip boundary conditions are the boundary conditions that we have discussed in the context of slip in this particular course. We have not discussed higher order slip boundary conditions 0.1 to 10 is typically the transition regime where it is a transition of the behavior from the continuum to the free molecular regime and typically one can capture the physics by using an equation which is a bit higher order accurate than the Navier-Stokes which is called as Bernat equation and with higher order slip boundary conditions. However you can also use the direct simulation Monte Carlo and the lattice Boltzmann for Knudsen number greater than 10 it is a free molecular flow. So one can use collision less Boltzmann equation DSMC or lattice Boltzmann model. Now we will consider molecular dynamics which is the agenda for this particular lecture. Molecular dynamics is essentially based on molecular behavior on the basis of Newtonian mechanics. So concept wise it is very very simple. So the idea if you say that so far as fluid mechanics goes I mean one has to carefully understand that molecular dynamics can also deal with quantum theories. But in the fluid dynamics context at least the quantum theories are not very relevant and we are interested about basically classical mechanics based molecular dynamics. So there are molecular dynamics simulations typically people working with chemistry often work with that kind of molecular dynamics simulations people working with quantum chemistry for example. So to get some properties actually one may start with quantum mechanics level calculations but for fluid mechanics we usually will not require that. So our basic law that governs the molecular behavior is the Newton's laws of are the Newton's laws of motion. So the idea if you say who introduce the idea of molecular dynamics it is nothing but Newton but obviously did not give the name molecular dynamics and its implementation however has been difficult because if see this is like a classical many body problem. So if there are many interacting bodies in fact more than 2 interacting bodies it is very difficult to solve the problem computationally. So in the computing edge with advancement in computational powers actually this is how like this is very interesting as I have showed that with advancement in manufacturing technology with advancement in fabrication technology microfluidics and nanofluidics has advanced. Similarly with advancement in computational power computational microfluidics and nanofluidics has also advanced. So the kinds of molecular dynamics simulations which were inaccessible even 10 years back with the present day highly powerful supercomputers it is I mean those calculations have become possible. So not that people did not know the principles but there were not enough platforms to implement that. So implementation has only been possible in the computing edge. These are computationally intensive the advantage of molecular dynamics is its theoretical background is very simple. Theoretical background I will discuss about this as simple as you can think just like elementary physics if you know and elementary numerical integration you know that is good enough. Post processing may require some statistical interpretation but solving the problem actually I mean its conceptually very simple but disadvantage is that you have limited system size. So you cannot simulate very large number of atoms or molecules. From atom positions velocities and accelerations calculate molecular positions and velocities at the next time step and integrate this infinitesimal steps to yield the trajectory of the system. And then there are efficient methods for integrating these elementary steps with several algorithms which we will discuss. When is molecular dynamics simulation most suitable? In considering MD simulation we need to keep the following things in mind. Simulation size. Simulation size by molecular dynamics we mean it is not a grid size it is the number of particles that we think of to extract observables over continuum length scales and also to minimize the boundary effects we need to consider many molecules. However limited computational resources mean that we can use at most a few hundred to thousand molecules I mean in normal computational platforms I am not talking about very highly powerful systems but even like I am talking about the systems which people commonly have access to even in laboratory environments which are not very rich. So, there you can get you you can do this. Thus it is inefficient for gas flows where large domains are needed due to large intermolecular distances and it deals most effectively with nano domains. The total time duration must be selected so that the calculation can finish with a reasonable time period but the simulation must be long enough to be relevant to the time scales of the processes being studied. Similarly we consider time scales of a few nanoseconds to microseconds so you can see that there are many physical processes which have time scales much larger than this. So molecular dynamics still is not good enough to capture the essential physics of those processes. The time scale the time step needs to be small enough to minimize discretization errors and also large enough for computation to be feasible. Now try to start with a very simple paradigm let us say you have a pressure driven flow in a nano channel. So if you have a pressure driven flow in a nano channel the continuum predictions will tell that this will be the velocity profile the fully developed velocity profile. However, with reduction in scale the exchange of momentum between the interfacial fluid layers and wall is incomplete leading to slip at the walls. This slip is often empirically accommodated in the form of a slip length. However you can directly obtain the velocity profile and hence the slip length from the molecular information from a force balance that on the molecular scale that gives a more realistic about idea about the actual flow mechanics in nano channels. So essentially that will give you a slip length which you can fit with continuum level calculations. For solving the molecular dynamics problem a particle based analysis is done by tracking the position and velocities of each molecule in the system. Now in the molecular scale simulations the deterministic path is the molecular dynamics where you solve equations of motion and integrate them numerically using small time steps. However in statistically based molecular simulations you can use probability based collision techniques combined with random movement of particles at successive time steps and one such technique well known technique is Monte Carlo simulation technique which we will not discuss in this particular course. So some basics of MDS. So you can see here that we can start with the Newton's second law of motion for each particle. When we say particle it is quote unquote particle so it could be atom molecule whatever. So particle does not necessarily mean particle in the context in the meaning of particle mechanics. So it could be atom molecule whatever. So now when you write the resultant force I will show you that we actually do not write the force directly. We write it as a negative of gradient of some potential because potential gives the interaction between molecule with another molecule or an atom with another atom much more explicitly. There are attractive potentials there are repulsive potentials and these can be directly captured by the potential rather than the force. But of course you can get the force from the potential. So then you can get the velocity and from the velocity so one integration of this once will give the velocity I will show the integration schemes and then you integrate it once more you get the position vector. So for a single particle you can get classical you can get analytical solution of this. In fact that is what we commonly do in high school physics. So you have a single particle you have various forces acting on the particle you find the position of the particle as a function of time. For n particle system it is a classical many body problem and with some empirical force field this problem cannot be analytically solved and that is where we use the computer power to solve this numerically. So in place of force we tend to use a gradient of a potential and the potential is always an approximation but it in physics what is the essential feature of modeling is we want to give a particular potential which represents or which mimics the physical reality to a good extent. In fact one important aspect of molecular modeling or molecular simulation is to design the potential because nowadays you have different software tools. You have different open source software which are which you can use and they are very reliable and which you can use for molecular dynamic calculation. So writing a code for molecular dynamics may not be a big issue but finding identifying a physically meaningful potential that capture the physics of the problem that you are considering is a matter of great interest. So designing potentials beyond the standard potential forms can be a matter of great interest in molecular modeling and molecular simulation not just using standard potentials but designing potentials. So I will first give a summary of the molecular dynamics and then we will discuss about the individual steps in details. So first we initialize the system and ensure that particles do not overlap in initial positions and we can randomly assign velocities. I will show you what how do you randomly assign velocities later on but this is just a summary. Then compute the energy. So if you have a system where chemistry is dominating then you can have an energy due to bond torsion, bond angle and non bonded form of energy which is beyond the chemical interaction. So I will show you different forms of interaction potential. Now you can write f is equal to minus dE dr is equal to m into a. Now you can write velocity from acceleration by noting this just a simple finite difference v at t plus delta t by 2 is equal to v at t minus delta t by 2 plus a into delta t. So it is just like a central difference type of scheme. So you can update v at t plus delta t by 2. So you have merged from time t to t plus delta t by 2 and then so half the time step you have merged and then the remaining half you merge as you find out r at t plus delta t is equal to r at t plus v at t plus delta t by 2 into delta t. So you have got the velocity at the midpoint of the time interval and using that midpoint you have as if jumped by another delta t by 2 and you have got the r at t plus delta t. So this algorithm is called as leapfrog algorithm very commonly used in molecular dynamics of very simple. I mean one of the very elementary steps infinite difference that you can use for this. The force fields where the physics comes or the chemistry comes. Many physical experiments were performed with gelatin or metal balls representing interacting liquid molecules. First computer experiments were performed using a hard sphere model but particles only interact during collisions and respond by an infinite repulsive force with no attractive component. So this is the potential energy positive potential energy will mean repulsive and the negative interaction energy will mean at attractive force. So you can see there is a sharp peak as you have hard spheres colliding soft sphere model but in reality you know atoms are not hard spheres. So soft sphere models improved over these models by replacing the infinite repulsion with a smoothly increasing repulsion and the pair potentials which were introduced later on try to most realistically represent the van der Waals interactions by combining the repulsive and the attractive components. So to compute the attractive and the repulsive components you can see that this the hard sphere model is the straight vertical line. But a more practical model where you have that as you find out the potential as a function of distance you have first the repulsive potential it comes down come to a minima and then like as you go beyond the critical distance the potential starts becoming attractive potential. And in the limit as you go to infinite r this attractive potential goes to 0. In reality you do not have to go for infinite r because calculation with infinite distance is computationally very very involved. So there is a cutoff distance typically a cutoff distance like you can see typically like this distance beyond which the potential goes to 0. So for gases you have no long range interactions and they interact by elastic collisions. So that means the hard sphere potential is not a bad choice for dense gases. So this is mainly for non-dense systems for dense gases liquids and solids you have electric dipole interaction which scales with r to the power – 3 you have if you have 2 dipoles that will scale with r to the power – 6. So the attractive component will be depending on that and you have repulsive nuclear forces which you may scale with r to the power – 12 I will show you what is the basis of r to the power – 12 and it is very very heuristic r to the power – 6 potential has a very sound mathematical platform and r to the power – 12 the platform is not so sound and I will show you what is the platform for that. So the total potential will have one term which is the attractive potential which will scale with r to the power – 6 and a repulsive potential which will scale with r to the power 12 and these 2 together is known as Lennard-Jones 6-12 potential or the Lennard-Jones potential. So the Lennard-Jones potential is based on the consideration that a pair of neutral atoms or molecules is subject to 2 distinct forces in the limit of large separation and small separation. An attractive force at long ranges which is typically Van der Waals force or dispersion force we have discussed about this earlier and a repulsive force at short ranges which comes into the picture as a result of overlapping electron orbitals referred to as Pauli repulsion from Pauli's exclusion principle. The Lennard-Jones potential also referred to as the LJ potential 6-12 potential or less commonly 12-6 potential is a simple mathematical model that represents this behaviour it was proposed in 1924 by John Lennard-Jones. So this is the mathematical form of the Lennard-Jones potential we have discussed about this here epsilon is the depth of the potential well sigma is the distance at which the inter particle potential is 0 and r is the distance between the particles. These parameters can be fitted to produce experimental data or deduce from results of accurate quantum chemistry calculations. So this is a very important thing where will you get sigma and epsilon. So to get these parameters actually one has to do quantum chemistry calculations. However one can also have controlled experimental data and one can fit that with this. The positive term describes repulsion and the negative term describes attraction which we have already discussed. The Lennard-Jones potential is an approximation. The form of the repulsion term has no theoretical justification. The repulsion force should depend exponentially on the distance but the repulsion term of the Lennard-Jones formula is more convenient due to the ease and efficiency of computing r to the power 12 as square of r to the power 6. This is just the basis. I mean it might appear to be very funny and unscientific but actually it works. So its physical origin is related to Pauli's principle. So that we have already discussed. Now Lennard-Jones attractive potential of physical basis that is the r to the power 6 component. Instantaneous position of electrons about nuclear protons may give rise to instantaneous dipole moments even if the atom is non-polar. We have discussed about this in dipole-dipole interaction the Van der Waals interaction. The instantaneous dipole generates an electric field that polarizes any nearby neutral atom inducing a dipole moment on it. The resulting interaction between two dipoles give rise to instantaneous attraction between two atoms with a finite time average. So as an example consider a Bohr atom in which one electron orbits around one proton. The first Bohr radius is defined such that E square bar by 4 pi epsilon naught A naught is equal to 2 H nu which gives A naught is equal to 0.053 nanometer. Here H nu is the energy of an electron in the first Bohr radius which is same as the energy needed to ionize the atom or the first ionization potential. The instantaneous dipole moment is A0 into E. Now consider a fixed dipole with a dipole moment mu interacting with an induced dipole of dipole moment mu induced. Strength of mu induced is dictated by the polarizability and is related to the electric field as the induced dipole moment is equal to alpha times E where alpha is the polarizability it is a constitutive equation. And the dipole moment basically if you have instantaneous charges of plus E and minus E then it is L into E where L is the perpendicular distance between the two. And expression for alpha for the Bohr atom may be derived by equating the external force on an electron with an internal force due to columbic interaction with the proton resolved along the direction of E. So you can see that you have E square sin phi where phi is this angle. So basically this interaction between plus E and minus E gives rise to a columbic force E square by 4 pi epsilon not S square that resolved in the direction of E is that into sin theta or sin phi. Now you can write sin phi as L by A0 so you can come up with this expression and this is equal to E into because E small e into L is mu induced. So this you can write E into mu induced by 4 pi epsilon not A0 cube okay. Now what is the external force? The external force is small e into capital E right. So f external is equal to f internal and that gives alpha is equal to 4 pi epsilon not A0 cube. Now what is E? If you have a dipole moment of mu this is an expression for E in terms of mu. This expression is not so important for us. What is important for us is that the details of the expression is not important but E scales with r to the power minus 3. So if E scales with r to the power minus 3 then the total interaction potential is minus half alpha E square. So E square scales with r to the power minus 6 and averaging over the angle theta when the angle theta is gone after integration you have basically V scales with E square and E scales with r to the power minus 3. So V scales with r to the power minus 6. This is the basis of the r to the power minus 6 in the Lennard Jones potential formula. This is very important because many times people use this as a blind formula but one should not do that and at least try to appeal to the physics which leads to this formula. Truncated Lennard Jones potential to save the computational time Lennard Jones potential is often truncated beyond r equal to rc is equal to 2.5 sigma this is a typical cut off radius. The corresponding Lennard Jones potential is about one sixtieth of its minimum value. So you can set the computational Lennard implementation of the Lennard Jones potential by using this formula for r less than rc and for r greater than rc. Now in the potential model when one considers the chemical bonds so you can see that there are different actions which are possible which can model the chemical bonds. So this model of chemical bonds was proposed by Linus Pauling in the 1930s I mean and this is one of big landmarks in chemistry. Bond angles and lengths are almost the same and the energy model can be broken up into two parts one is the covalent terms and the non-covalent terms. Now what are the interactions that we are talking about the interaction forces. So you can have stretching of the bonds so bonds may act like springs so you can have stretching of the bonds you can have bending of the bonds and you can have torsion of the bonds. So torsion of the bond means basically twisting with respect to its axis that will give the torsion of the bonds and bending of the bonds means it is basically it is bending in the plane in which the bonds are present. So you can have the total energy as stretching energy the bonding energy and torsion energy and non-bonded forms of energy. So the stretching energy you can see that this stretching energy is typically like of the form k into r-r0 square this is like just like a spring with k0 is the spring stiffness. Model different bonds by playing with the parameter which is the stiffness. The bending energy is related to the bending angle theta. So it is very similar to the stretching energy but the parameters are different and the bending angles are more relevant here and not just the position vector okay this we will skip. The torsion energy you can see that you can represent the torsion energy also by a suitable mathematical form which can be say for example a periodic function by which you represent the torsion energy. For example in this formula the capital A controls the amplitude, N controls the periodicity and Phi controls shifts from the or shifts Phi shifts the entire curve along the rotation angle axis. So various parameters can control various interaction. So with different interaction parameters you can have different torsional energies. Now non-bonded form of energy. So you can have the van der Waals terms. So this is what is the attractive term and this is the repulsive term and then you have the electrostatic terms and electrostatic interactions are important in a system in which electrical field is there and ions are present. So electrostatic terms may be important and electrostatic terms you can see these are expressed just by using the Coulomb's law. So the non-bonded energy parameters the Lennart zones like just to visually give you an appeal that how the Lennart zones parameters the Lennart zones potential will vary if you vary the parameters. And these parameters are already important you know you can play with these parameters to mimic hydrophobicity and the weightability of the substrate. So these parameters if you alter the weightability of the substrate, different weightabilities of the substrate can be represented. So the physics of the weightability of the substrate can be translated to the molecular model by altering or tuning with the Lennart zones parameters, Coulombic interactions. So Coulombic interactions like you have this kind of force Q i Q j by 4 pi epsilon not R square this is the force due to Coulombic interaction. It decays as 1 by R square it is a long range force and its tail extends longer than the attractive Van der Waals force. Long range contributions need to be calculated by a technique called as Ewald summation. So Ewald summation of course there are other methods also but Ewald summation is the most commonly used technique. And what is done in Ewald summation the ith particle with charge Q i is surrounded by a diffuse charge distribution of opposite sign. The screen charge potential decays rapidly and thus becomes a short range force. However one needs to calculate the potential distribution for the screen charge. A compensating cloud for screening is added and the contribution is added to the screen charge potential. The compensating charge distribution is a Gaussian whose contribution has an analytically closed form. So point charges can be equivalently represented by point charges plus a screening cloud plus a compensating cloud that compensates the screening. Now in engineering we are many times interested in modeling water through nano channels. So how do you represent the interactions in water in a molecular dynamics framework. Intermolecular interaction of water molecules is expressed as a combination of Van der Waals force and electrostatic force and it is the electrostatic force is very very important because the water will have a dipolar behavior. So considering water as the fluid brings in the effect of hydrogen bonding that makes it so special. So modeling with water is actually not a very simple thing. Still we do not understand water very well and in molecular dynamics you will see there are many papers even papers very sort of highly contributing papers are these days published with reduced substances with substances with very special properties like properties of inert gases like that inert substances. So not many studies are reported where you have dynamics include flow dynamics which addresses water because water is not so simple to deal with in a non equilibrium environment. So you can have HOH bond angle for water typically 109.47 degree that is taken in a model we call which we call as SPCE model simple point charge extended model which includes the effect of permanent and induced dipoles considering the polarization effects. So in the SPCE model you have Leonard Jones standard Leonard Jones potential standard electrostatic potential plus a polarization potential. So these together make the SPCE model which is the possibly the most commonly used models for potential of water. So overall kinetics you start with the force field which is negative of the potential field you update the so you can update the position and velocity by using the force field and you can get the system temperature by relating half m v square is equal to 3 by 2 n k B T. So n half k B T for each degree of freedom so you will give total get as 3 by 2 n k B T. Sometimes molecular dynamics is typically some not sometimes very common very often molecular dynamics is typically applied to systems of a few 100 to 1000 atoms such small systems are usually dominated by surface effects. If we want to simulate the bulk liquid the surface effects are removed using periodic boundary conditions if n molecules are confined to a volume v we imagine that it is surrounded by exact replicas of itself. We do not need to follow trajectories of the image molecules because they can be easily computed when needed each molecule in the primary cell interacts with all the n-1 other molecules. So the molecule moves into the image cell the image from opposite cell moves into the primary cell so it is basically repetitive behavior of some monoblocks where each monoblock is a unit cell type of thing and you basically consider that whatever is the boundary condition at the inlet for example at the outlet of the repeating unit the same boundary condition is repeated. So that is called as periodic boundary condition initialization. So for implementation this is what is important it involves assignment of initial position velocities and higher order derivatives of positions. In equilibrium molecular dynamics the independent thermodynamic properties are n, v and e typically n is between 100 to 1000 so we can capture the phenomena of interest and keep the simulation feasible. Given n, v and e we can compute the number density of molecules and the energy per molecule. The initial velocities are randomly generated with the constraint that the total kinetic energy is related to the given temperature and we want the simulation to have or we can choose the total energy and set the positions and velocities to correspond. The initial distribution of velocities is determined from a random distribution with the magnitudes conforming to the required temperature and corrected so that there is no overall momentum. So this is a very important constraint it needs to be satisfied the velocities are chosen randomly from an equilibrium distribution known as Maxwell Boltzmann distribution which gives the probability that a molecule i has a velocity vx in the x direction at a temperature t. The positions may be taken also from a previous simulation. In the case of initial positions are created at random this must be relocated with energy minimization technique with minimum potential energy of the system so as to render the system realistic. So then over the first few hundred steps the system relaxes from the initial conditions and approaches equilibrium. The relaxation phase is called as equilibration. It consists of three steps apply finite difference algorithm accumulate monitored properties scale the velocities to keep it at the required temperature. So if you are basically simulating an isothermal system the duration of the equilibrium varies depending on how far remove the initial conditions are from the equilibrium state. So initially your initial configuration should correspond to an equilibrium state. So that equilibrium state if your initial choice is far from equilibrium then several steps will be required for equilibration. How do we know when equilibrium is attained? This is a very important issue we know from kinetic theory that the velocity distribution at equilibrium should be Maxwell's Gaussian distribution. As this is cumbersome to monitor the velocity distribution we monitor a single number Boltzmann's H function. So this is from statistical mechanics that this is a typical function where fvx is the actual velocity distribution defined as the fraction of molecules having velocity vx plus minus half delta vx. This is compared with the H function of the actual Maxwell's distribution as the simulation moves towards the equilibration the instantaneous H function should converge to the expected value. If the initial velocities are not random this will take longer. The numerical integration I have already discussed about this so there are several techniques like Verlet technique velocity Verlet predictor character several other techniques. So like the velocity rt plus delta t plus rt minus delta t by 2 delta t into r double dot t sorry this is the average velocity. So r at t plus delta t plus r at t minus delta t by 2 delta t so this you basically get by using a finite difference type of formula. So it is like it is basically so there is a plus minus sign problem here I think. So basically the difference r at t plus delta t minus r at t minus delta t so not plus divided by 2 delta t. So it is just like a finite difference calculation then velocity Verlet the velocity calculated at mid step and positions for velocities are improvements over the Verlet with a greater accuracy. So this is just like the leapfrog type of algorithm that I have discussed earlier. Next concept is a very important concept if you are implementing an isothermal problem thermostats. For several reasons example drift during equilibration drift as a result of force truncation and integration errors heating due to external or frictional forces it is necessary to control the temperature of the system in molecular dynamic simulations. So you have different algorithms which brings the temperature back to the original temperature which has artificially shifted from its value for a problem which is actually an isothermal problem. For example I am talking about one algorithm which is called as Berenson algorithm what it does is it mimics weak coupling with first order kinetics to an external heat bath with a given temperature T0. The effect of this algorithm is that a deviation of the system temperature from T0 is slowly corrected. So you employ a dt dt where small t is the time is equal to T0 minus t by tau where tau is the time constant. So this means that the temperature deviation decays exponentially if you integrate this then the temperature versus time will be exponentially decaying. So temperature deviation decays exponentially with a time constant tau okay. So there are very involved algorithms for many thermostats but essentially the thermostat what it is trying to do is that if there is a deviation from the design temperature it attempts to bring back the temperature to the design temperature over a period of time through an exponential decay of the temperature gap. Now you get molecular information in molecular dynamics but how do you determine macroscopic properties? So here you use some statistical mechanics. So it has to be to determine a macroscopic property from phase space trajectories. So the phase space will typically have the position the velocity versus position or momentum versus position space. So you can see sorry so you can see that it is a function of the position and momentum or you can write it position or velocity position and velocity. So if delta t is the time step you basically use the average or expectation value to get the property A. This is an approximation of the infinite time average. So what you do is that instead of using this averaging that is integration of the property over the time you simply do a summation. You replace the integration with a summation in the phase space that is what you do. So the overall flow chart you have an initial configuration and potential based on that you calculate force, obtain velocity and then obtain position and using statistical mechanics you do a post processing. So when you do the post processing and then if you want to extract the physical units see in molecular dynamics we commonly deal with the Lennard-Jones units. So one unit of length, one unit of temperature all these things in the molecular dynamics paradigm is different from physical units. So how do they correspond to physical units? So the reduced units are commonly used for simulating Lennard-Jones fluids. You set sigma equal to epsilon is equal to m equal to kBt equal to 1. Now the length is normalized by the parameter sigma mass by the molecular mass time by this parameter. So you can see these are normalized like r into sigma to the power minus 1 this is a normalized unit. So similarly the time, temperature, energy I will give you one example of pressure like the pressure is P into sigma cube epsilon to the power minus 1. So one unit of pressure has to be multiplied with epsilon by sigma cube to convert to the corresponding physical unit of pressure. So one applying one unit of pressure in a Lennard-Jones unit in a sample of water molecules will mean that the actual pressure physical pressure is 1 into epsilon by sigma cube. So then if you use that is the typical value of epsilon then that is 0.65 kJ per mole this is the basically the interaction energy parameter and sigma 0.3166 nanometer I mean typical length scale of the water molecule then that is roughly 10 to the power 8 Pascal. So many times in molecular dynamics we say one unit of pressure that will typically mean this kind of physical unit. So as an engineer we should be comfortable with converting from the reduced unit to the physical unit because we are commonly asked to design a practical problem where the units are physical units are not just the Lennard-Jones or the reduced units. So you can so post process the properties that you get from the molecular dynamics simulation. So some of the static properties and the dynamic properties so you can post process a wide range of properties I am not going into all this. So what are the concepts that we assimilated from the molecular dynamics simulation? The molecular flows and the tools, molecular dynamics simulations the inputs, outputs, force fields and the algorithms and the post processing. Now I would like to conclude this particular lecture by discussing about where is the future direction of research going in these areas. Of course there are many physical problems we are trying to capture the physics of the problem in the nanoscale physics of flow in the nanoscale through molecular dynamics. This is a pure sort of physical understanding of the scenario in the molecular scale. Computational issues, there are computational issues associated with the possibility of having massively parallelized systems where you can implement a large number of molecules and so on that is one side. There is another very interesting side which is called as multi-scale modeling. So let us say that you have a nano channel where you want to solve the problem by using molecular dynamics simulations. Even if it is nano channel it may actually contain a large number of molecules at least large than few hundreds of thousands. So it is very difficult to compute the system with that many number of molecules. So what you essentially can do is close to the wall where molecular interactions are important you go for molecular dynamics simulations with a few thousands of molecules. In the bulk you use still the Navier-Stokes equation maybe with certain modifications if necessary and then you patch up these two. So there may be overlapping layer where you have the domain of the molecular dynamics overlapping with the domain of the continuum simulation and they are exchanging information with each other. So in this way in the same simulation paradigm you can capture multiple physical scales. This is called as multi-scale modeling. Many of you are interested in CFD and I can tell you that in the classical sense of CFD there is very little new aspect of research that one can undertake at this moment. The main new aspects of research in the classical CFD are based on solving newer and newer physical problems but not in the algorithm itself. Now where you have a big challenge in CFD is this multi-scale CFD and multi-scale CFD with multiple physics connected with it. For example physics of electricity and magnetism connected with fluid mechanics. So then the modern direction of CFD where it is going towards is a combination of multi-scale and multi-physics CFD. So a CFD that can capture multiple physical issues over multiple physical scales and molecular dynamic simulation in conjunction with the continuum based calculations can take one effectively towards that direction. A good amount of research is being done currently in this area but still there are many open ended issues which remain to be resolved in this particular area of interest. Thank you very much. We stop here today and in the next lecture we will start with a new topic.