1 / 26

Computational Physics (Lecture 5)

Computational Physics (Lecture 5) . PHY4370. Extreme of a function. An associated problem to finding the root of an equation is finding the maxima and/or minima of a function . Examples: the equilibrium position of an object, the potential surface of a field,

medea
Download Presentation

Computational Physics (Lecture 5)

An Image/Link below is provided (as is) to download presentation Download Policy: Content on the Website is provided to you AS IS for your information and personal use and may not be sold / licensed / shared on other websites without getting consent from its author. Content is provided to you AS IS for your information and personal use only. Download presentation by click this link. While downloading, if for some reason you are not able to download a presentation, the publisher may have deleted the file from their server. During download, if you can't get a presentation, the file might be deleted by the publisher.

E N D

Presentation Transcript

  1. Computational Physics(Lecture 5) PHY4370

  2. Extreme of a function • An associated problem to finding the root of an equation is finding the maxima and/or minima of a function. • Examples: the equilibrium position of an object, the potential surface of a field, and the optimized structures of molecules small clusters and crystals.

  3. Condition of extreme • Condition: f (x) = dg(x)/dx= 0, • f’ is greater (less) than 0 => minimum (maximum) • Therefore, searching of root can be generalized here to search for the extremes. • At each step of updating the value of x, • make a judgment as to whether g(xk+1) is increasing (decreasing) • if we are searching for a maximum (minimum). • If it is, we accept the update. • If it is not, we reverse the update; that is, instead of using xk+1 = xk+ Δxk, we use xk+1= xk− Δxk. • Newton method, the increment is Δxk= −fk/ fk’ • the secant method, the increment is Δxk= −(xk− xk-1) fk/( fk− fk-1).

  4. One simple example • finding the bond length of molecular NaClfrom the interaction potential between the two ions • Assuming that the interaction potential is V(r ) when the two ions are separated by a distance r , the bond length req is the equilibrium distance when V(r ) is at its minimum. • The interaction potential is: • V(r ) = − e^2/4πε0r+ V0exp(−r/r0) • where e is the charge of a proton, ε0is the electric permittivity of vacuum, and V0 = 1.09x10^3 eV and r0 = 0.330 A are parameters of this effective interaction. • The first term: Coulomb interaction, • The second term: result of the electron distribution in the system.

  5. At equilibrium, the force is zero • Simply search for the root of f(x) = d g(x) /dx=0 (g(x) = -V(x)) • Sample code: use secant method to find the bond length of NaCl.

  6. // An example of calculating the bond length of NaCl. import java.lang.*; public class Bond { static final double e2 = 14.4, v0 = 1090, r0 = 0.33; public static void main(String argv[]) { double del = 1e-6, r = 2, dr = 0.1; int n = 20; r = secant2(n, del, r, dr); System.out.println("The bond length is " + r + " angstroms"); } // Method to carry out the secant search for the // maximum of g(x) via the root of f(x)=dg(x)/dx=0. public static double secant2(int n, double del, double x, double dx) { int k = 0; double x1 = x+dx, x2 = 0; double g0 = g(x); double g1 = g(x1); if (g1 > g0) x1 = x-dx; while ((Math.abs(dx)>del) && (k<n)) { double d = f(x1)-f(x); dx = -(x1-x)*f(x1)/d; x2 = x1+dx; double g2 = g(x2); if (g2 > g1) x2 = x1-dx; x = x1; x1 = x2; g1 = g2; k++; } if (k==n) System.out.println("Convergence not" + " found after " + n + " iterations"); return x1; } // Method to provide function g(x)=-e2/x+v0*exp(-x/r0). public static double g(double x) { return -e2/x+v0*Math.exp(-x/r0); } // Method to provide function f(x)=-dg(x)/dx. public static double f(double x) { return -e2/(x*x)+v0*Math.exp(-x/r0)/r0; } }

  7. The bond length obtained from the above program is 2.36 A. • One important issue is that in many cases we do not have the explicit function f (x) = g’(x) • if g(x) is not an explicit function of x. • We can construct the first-order and second-order derivatives numerically.

  8. Steepest decent method • In the sample program above, the search process is forced to move along the direction of descending the function g(x) when looking for a minimum. • xk+1= xk+ Δxk, the increment Δxkhas the sign opposite to g’(xk). • An update scheme can be formulated as • Xk+1= xk+ Δxk= xk− ag’(xk), • a is a positive, small, and adjustable parameter. • In a multivariable case • Xk+1= xk+ Δxk= xk− a∇g(xk)/|∇g(xk)|, • where x = (x1, x2, . . . , xl) and ∇g(x) = (∂g/∂x1, ∂g/∂x2, . . . , ∂g/∂xl ). • Note that step xkhere is scaled by |∇g(xk)| and is forced to move toward the direction of the steepest descent. => steepest descentmethod. • One example: an implementation of such a scheme to • search for the minimum of the function g(x, y) = (x − 1)^2 exp(−y^2) + y(y + 2) exp(−2x^2)around x = 0.1 and y = −1.

  9. // An example of searching for a minimum of a multivariable // function through the steepest-descent method. import java.lang.*; public class Minimum { public static void main(String argv[]) { double del = 1e-6, a = 0.1; double x[] = new double[2]; x[0] = 0.1; x[1] = -1; steepestDescent(x, a, del); System.out.println("The minimum is at" + " x= " + x[0] +", y= " +x[1]); } // Method to carry out the steepest-descent search. public static void steepestDescent(double x[], double a, double del) { int n = x.length; double h = 1e-6; double g0 = g(x); double fi[] = new double[n]; fi = f(x, h); double dg = 0; for (int i=0; i<n; ++i) dg += fi[i]*fi[i]; dg = Math.sqrt(dg); double b = a/dg; while (dg > del) { for (int i=0; i<n; ++i) x[i] -= b*fi[i]; h /= 2; fi = f(x, h); dg = 0; for (int i=0; i<n; ++i) dg += fi[i]*fi[i]; dg = Math.sqrt(dg); b = a/dg; double g1 = g(x); if (g1 > g0) a /= 2; else g0 = g1; } }

  10. // Method to provide the gradient of g(x). public static double[] f(double x[], double h) { int n = x.length; double z[] = new double[n]; double y[] = (double[]) x.clone(); double g0 = g(x); for (int i=0; i<n; ++i) { y[i] += h; z[i] = (g(y)-g0)/h; } return z; } // Method to provide function g(x). public static double g(double x[]) { return (x[0]-1)*(x[0]-1)*Math.exp(-x[1]*x[1]) +x[1]*(x[1]+2)*Math.exp(-2*x[0]*x[0]); } }

  11. Note that the spacing in the two-point formula for the derivative is reduced by a factor of 2 after each search, so is the step size in case of overshooting. • After runningtheprogram, the minimum is at x =0.107 355 and y =−1.223 376. • simple but not very efficient scheme. • It converges to a local minimum near the starting point. • The search for a global minimum or maximum of a multivariable function is a nontrivial task, especially when the function contains a significant number of local minima or maxima. More will be covered later.

  12. An Overview on Inter-atomic Potentials • Simulation results are only as good as the potentials (forces) that are used. • It is not an easy problem to get accurate forces because the origin of forces is quantum mechanics, to date a computationally intractable problem, but there has been some recent progress. • The most common approach is semi-empirical: a functional form is guessed or deduced and parameters are fitted from experimental data or ab initio potentials calculated for a model system.

  13. The Born-Oppenheimer Approxiamtion • How do we define the potential to be used? • The fundamental description is in terms of the non-relativistic Schröedinger equation for both the ionic and electronic degrees of freedom. • But since the electrons are much lighter they respond more-or-less instantly to the ionic motion. • Also typical temperatures are very cold compared to the electronic energy scale (1 Hartree, which sets the scale of the electronic energy, equals 316,000 K) • so the electrons are in or close to their ground state.

  14. In the Born-Oppenheimer (BO) approximation we fix the ions (and neglect the nuclear kinetic energy) and solve (or try to solve) for the electronic ground state energy. • Corrections to BO are usually very small because corrections go as the ratio of the electron to nuclear mass which is smaller than 0.0001 for all atoms except hydrogen and helium.

  15. Unfortunately solving for the quantum mechanics is not very easy to do! • The most popular methods are the local density functional theory (and improvements such as the GGA theory), quantum chemistry methods such as Hartree-Fock theory and the configuration interaction method. Another (accurate) method is quantum Monte Carlo.

  16. One break-though occurred fifteen years ago when Car and Parrinello (CP) showed how one could simultaneously solve the LDA equation for the electronic energy and the dynamical equations for the nuclei. • This means one does not have to tabulate the BO energy but just compute it as the ions move. We will describe this method later in the course. But the CP approach is very much slower than using fitted potentials (about 4 orders of magnitude). The largest CP studies have several hundreds of electrons or atoms while the largest MD simulation have millions of particles.

  17. Unless we use the CP approach then we have to assume some form for the potential. • The assumption of a potential is a big approximation! • the potential energy surface is a function of 3N variables for N ions. • Symmetries such as translational degrees or rotational of freedom will only get rid of a few degrees of freedom. • What usually has to be done is to assume the many-body potential is a sum of pair potentials and perhaps three-body potentials etc. • This also assumes the molecule is rigid and weakly interacting with the environment so its internal electronic state does not couple with surrounding molecules. • A pair potential is a function of a single variable; there is no way it can exactly represent a 3N dimensional function.

  18. So where do we get the pair potentials from? • The basic input is theoretical as we will discuss. The other input is empirical. Some sources of experimental data are: • atom-atom scattering in the gas phase can give direct information about the atom-atom potential. The scattering energy samples different regions of the interaction depending on the energy of the colliding atoms. • gas phase data, for example the virial coefficients and transport properties give integrals over the potential. • the low temperature thermodynamic properties. Most materials (except helium) form a solid at low temperature. • If the energy as a function of density is known, that gives fairly direct information about the depth of the potential well and its curvature.

  19. high pressure information. • This samples the interaction at shorter distances than at lower pressures. • thermodynamic data, • such as melting temperature, the critical point and the triple point. This information is harder to use directly as it is many-body property. But one can use simulations to compute a phase diagram with a trial potential, then try to change the potential using perturbation theory to make the trial phase diagram agree with the actual phase diagram.

  20. The difference between interpolation and extrapolation. • If experimental information is used to determine the potential, it will be good in a nearby region of phase space and for similar properties. • But that potential could be quite bad for unrelated properties or outside the region that has been fitted. • For example, if the potential has been determined using bulk properties, there is no reason to think it will be good to compute surface properties because the local environment on a surface is not as symmetrical as in the bulk.

  21. Molecular Solids: The Noble Gases • These are the rare gases where the electronic shells are closed. • This means they are spherical atoms. They were the first systems to be simulated and the pair potentials work best for them! • At large distances one can understand the interactions by considering how two oscillators would interact, giving the r -6 attraction where the coefficient is determined by the atomic polarizability. • On the other hand at short distances, the potentials are strongly repulsive because of the Pauli exclusion principle (electrons are fermions). • Lennard-Jones assumed a stronger inverse power, conveniently r -12. One should not think the LJ potential is anything more fundamental than this.

  22. Even for Argon it is only accurate to 10% or so. And three-body potentials can be significant at the 10% level. • The potentials that are fit to bulk data are really effective potentials and could be somewhat different than the real 2-body potential that would be determined from gas-phase data because they include correction for three- and higher-body interactions. Hence the effective potentials can be dependent on density and less so on temperature.

  23. Metals • How about the metals? • Metals are characterized by a rare gas shell plus a number of valence electrons. • Pair potentials do not work terribly well because the electrons are very much delocalized. • If one does use a pair potential, it must be softer than the Lennard-Jones form, for example a Morse potential (sum of exponentials.) • The major problem with pair potentials for solid metals is that they have no environmental dependence so that a bulk atom is too much like a surface atom.

  24. One popular empirical potential for metals is the embedded-atom-method (EAM) potential. It is based on the basic idea of jellium model that “metallic bonding” involves ions in an electron sea (jellium) so the total potential energy is the sum of direct ion-ion interaction • (for example to keep to ion cores from overlapping) and the interaction between an ion and the electron sea. • One assumes that locally there is an electron charge density coming from the nearby ions and that the energy and interaction are some function of this density. • This later term is then a true many-body potential. • The EAM potential works well for spherically symmetric atoms such as Cu, Al, Pb. Some references are Daw and Baskes, Phys. Rev. B 29, 6443 (1984) andFoiles, Baskes and Daw, Phys. Rev. B 33, 7983 (1986). • Tight-binding potentials are another possible functional forms for treating metals, especially for transition-metal solids.

  25. Metal: EMA potential:

More Related