Mod-01 Lec-45 Introduction to Molecular Dynamics Simulations (Contd.)

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 a Navier stokes equation Like for example, one order more accurate than a Navier stokes equation, these are known as Burnett equations, 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 what 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 modelling 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 may be 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 these 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 that 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 base 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; < 10 to the power - 3 these are, again I; as I told these are not hard and first numbers, just typical orders you can use Navier stokes equation with no slip boundary condition at the wall between 10 to the power - 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 behaviour 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 Burnett 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 > 10, it is a free molecular flow, so one can use collisionless 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 behaviour 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 dynamic simulations Typically, people working with chemistry often work with that kind of molecular dynamic 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 behaviour is the Newton’s laws of or the Newton’s laws of motion So, the idea if you say, who introduced the idea of molecular dynamics, it is nothing but Newton but obviously, it 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 age 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, micro fluidics and nano fluidics has advanced. Similarly, with advancement in computational power, computational micro fluidics and nano fluidics has also advanced So, the kinds of molecular dynamic 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 age, these are computationally intensive The advantage of molecular dynamics; theoretical background is very simple; theoretical background I will discuss about is as simple as you can think just like elementary physics, if you know an elementary numerical integration you know that is good enough. Post processing may require some statistical interpretation but solving the problem actually, I mean it is 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 these 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, 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 100 to 1000 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 can do this Does 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 timescales of the processes being studied Usually, we consider timescales of a few nanoseconds to microseconds So, you can see that there are many physical processes, which have timescales much larger than these, 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, I will 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 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 these ones 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 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 softwares which are; which you can use and they are very reliable and which you can use for molecular dynamics calculations. So, writing a code for molecular dynamics may not be a big issue but finding, identifying a physically meaningful potential that captures 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, it is 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 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 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 = – dE/

dr is = m * a, now you can write velocity from acceleration by noting this just a simple fine difference, v at t + delta t/ 2 is = v at t – delta t/ 2 + a * delta t, okay, so it is just like a central difference type of scheme. So, you can update v at t + delta t/2, so you have merger from time t to t + delta t/ 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 + delta t is = t at t + v at t + delta t/ 2 * 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/ 2 and you have got the r at t + delta t, so this algorithm is called as leapfrog algorithm, very commonly used in molecular dynamics. So, very simple, I mean one of the very elementary steps in finite difference that you can use for this The force fields, where the physics comes or the chemistry comes historically, physical experiments were performed with gelatin or metal balls representing interacting liquid molecules. First computer experiments were performed using a hard sphere model, particles only interact during collisions and respond via 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 attractive force, so you can see there is a sharp peak as you have hard spheres colliding, soft sphere; but in reality you know atoms are not hard spheres, so software 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 at 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 a 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 a very very involved, so there is a cut off distance typically, a cut off 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 this two 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 And 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 exclusion principle. The Lennard Jones potential also referred to as the LG 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 2l 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, a 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 nonpolar, 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 2 dipoles give rise to instantaneous attraction between 2 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/ 4 pi epsilon 0 a0 is = 2h Nu, which gives a0 = 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 * e Now, consider a fixed dipole with a dipole moment Mu interacting with an induced dipole of dipole moment Mu induced, okay. 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 +e and –e, then it is l * e, where l is the perpendicular distance between the 2 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 + e and – e give rise to a columbic

force e square/ 4 pi epsilon 0 a square that resolved in the direction of capital E is that into sin theta or sin phi Now, you can write sin phi as l/a0, so you can come up with this expression and this is = e *; because small e * l is Mu induced, so this you can write e * Mu induced/ 4 pi epsilon 0 a0 cube okay. Now, what is the external force? The external force is small e * capital E, right, so F external is = F internal and that gives alpha is = 4 pi epsilon 0 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 – 3, so if E scales with the r to the power – 3 then the total interaction potential is – 1/2 alpha e square. So, E square scales with r to the power – 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 – 3, so V scales with r to the power – 6, this is the basis of the r to the power – 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 = rc is = 2.5 sigma, this is a typical cut off radius, the corresponding Lennard Jones potential is about 1/ 60 th of its minimum value. So, you can set the computational Lennard; implementation of the Lennard Jones potential by using this formula for r < rc and for r > 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 1930’s 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 2 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, a stretching energy, the bonding energy and torsional 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 * r – r0 square this is like just like a spring with k0 is the spring stiffness, so you can 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 a 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 that torsion energy. For example, in this formula that capital A controls the amplitude, n control the periodicity and phi control shifts from the; or shifts; phi shifts the entire curve along the rotation angle axis So, various parameters can control various interactions, so with different interaction parameters, you can have different torsional energies Now, non bonded form of energy; so you can have the van der Waal’s 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 Lennard Jones like just to visually give you an appeal that how the Lennard Jones parameters; the Lennard Jones potential will vary, if you vary the parameters and these parameters are very important you know, you can play with these parameters to mimic hydrophobicity and the weight ability of the substrate, so these parameters if you alter, the weight ability of the substrate; different weight abilities of the substrate can be represented So, the physics of the weight ability of the substrate can be translated to the molecular model by altering or tuning with the Lennard Jones parameters. Coulombic interactions; so coulombic interactions like you have this kind of force qi, qj/ 4 pi epsilon 0 r square, this is the force due to coulombic interaction It decays as 1/ r square, it is a long range force and its tail extends longer than the attractive van der Waal’s force Long range contributions need to be calculated by a technique called as a 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 qi is surrounded by a diffuse charge distribution of opposite sign. The screened 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 screened 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 Vander Waals force and electrostatic force and it is the electrostatic force is very very important because the water will have dipolar behaviour 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 what are 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; flow dynamics which addresses water because water is not 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 Lennard Jones; standard Lennard 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 1/2 mv square is = 3/2 N kv T, so 1/2 kv T for each degree of freedom, so you will give total gate as c/ 2 N kv T Sometimes, molecular dynamics is typically; 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. If a molecule moves into the image cell, the image from opposite cell moves into the primary cell, so it is basically repetitive behaviour of some mono blocks, where each mono block 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 Initialisation; 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 these 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 100 steps, the system relaxes from the initial conditions and approaches equilibrium. The relaxation phase is called as equilibrium, it consists of 3 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 removed the initial conditions are from the equilibration state; 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 f vx is the actual velocity distribution defined as the fraction of molecules having velocity vx +- 1/2 delta vx This is compared with the H function of the actual Maxwell distribution, as the simulation moves towards 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 corrector, several other techniques, so like the velocity rt + delta t + rt – delta t /2 delta t * r double dot t; sorry this is the average velocity, so r at t + delta t + r at t – delta t/ 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, it is basically the difference r at t + delta t – r at t – 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 equilibrium 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 a different algorithm, 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 Berendsen 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 = T0 – T/ tau, where tau is the time constant So, this means that the temperature deviation decays exponentially, if you integrate this then that 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 = epsilon is = m = k BT = 1, okay, 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 * sigma to the power -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 * sigma cube epsilon to the power – 1 So, one unit of pressure has to be multiplied with epsilon/ sigma cube to convert to the corresponding physical unit of pressure. So, one applying one unit of pressure in the Lennard Jones unit in a sample of water molecules will mean that the actual pressure physical pressure is 1 * epsilon/ sigma cube, so then if you use that is the typical value of epsilon then that is 0.65 kilo joule per mole, this is the basically the interaction energy parameter And sigma 0.3166 nano meter, 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 and 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 these 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 nano scale; physics of flow in the nano scale 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 parallelised 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 modelling, 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 100’s or 1000’s, 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 1000’s of molecules. In the bulk you use still a Navier Stokes equation maybe with certain modifications, if necessary and then you patch up these 2, 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 modelling. 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 dynamics 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