CE Molecular Simulation

Document Sample
CE Molecular Simulation Powered By Docstoc
					                                              1




CE 530 Molecular Simulation

               Lecture 24
  Non-Equilibrium Molecular Dynamics

                   David A. Kofke
         Department of Chemical Engineering
                   SUNY Buffalo
               kofke@eng.buffalo.edu
                                                                  2




                             Summary
                            from Lecture 12
 Dynamical properties describe the way collective behaviors
  cause macroscopic observables to redistribute or decay
 Evaluation of transport coefficients requires non-equilibrium
  condition
   • NEMD imposes macroscopic non-equilibrium steady state
   • EMD approach uses natural fluctuations from equilibrium
 Two formulations to connect macroscopic to microscopic
   • Einstein relation describes long-time asymptotic behavior
   • Green-Kubo relation connects to time correlation function
 Several approaches to evaluation of correlation functions
   • direct: simple but inefficient
   • Fourier transform: less simple, more efficient
   • coarse graining: least simple, most efficient, approximate
                                                                       3




    Limitations of Equilibrium Methods
 Response to naturally occurring (small) fluctuations
 Signal-to-noise particularly bad at long times
   • but may have significant contributions to transport coefficient
     here
 Finite system size limits time that correlations can be
  calculated reliably




          correlations between                  …lose meaning once
          these two…                            they’ve traveled the
                                                length of the system
                                                                 4




  Non-Equilibrium Molecular Dynamics
 Introduce much larger fluctuation artificially
   • dramatically improve signal-to-noise of response
 Measure steady-state response
 Corresponds more closely to experimental procedure
   • create flow of momentum, energy, mass, etc. to measure…
   • …shear viscosity, thermal conductivity, diffusivity, etc.
 Advantages
   • better quality of measurement
   • can also examine nonlinear response
 Disadvantages
   • limited to one transport process at a time
   • may need to extrapolate to linear response
                                                                         5




            One (Disfavored) Approach

 Introduce boundaries in which molecules interact with
  inhomogeneous momentum/mass/energy reservoirs
 Disadvantages
   •   incompatible with PBC
   •   introduces surface effects
   •   inhomogeneous
   •   difficult to analyze to obtain transport coefficients correctly
 Have a look with a thermal conductivity applet
 Better methods rely on linear response theory
                                                                                           6




       Linear Response Theory: Static

 Linear Response Theory forms the theoretical basis for
  evaluation of transport properties by molecular simulation
 Consider first a static linear response
 Examine how average of a mechanical property A changes in
  the presence of an external perturbation f
   • Unperturbed value A 0
   • Apply perturbation to Hamiltonian H  H0   B( p N , q N )
                                    d Ae  0
                                           H  B 
   • New value of A               
                         A 0  A                  H0  B 
                                         d e                     Susceptibility
   • Linearize
                   (A)                                         describes first-order
                              AB
                                            A       B 0
                                                         
                     0            0         0                static response to
                                                                   perturbation
                                                                           7




     Example of Static Linear Response

 Dielectric response to an external electric field
   • coupling to dipole moment of system, Mx       H   E y M y (q N )
   • interest in net polarization induced by field M y
   • thus A = B = My
            No field                          Field on
                                                              E



            My  0                             My  0
                       2                2
   • susceptibility   M y
                             0
                                   My
                                         0
                                          
                                                                                              8




   Linear Response Theory: Dynamic 1.

 Time-dependent perturbation Fe(t)
 Consider situation in which Fe is non-zero for t < 0, then is
  switched off at t = 0
 Response A decays to zero
                                       H 0  B 
           A(t )     d A(t )e                      Ensemble average over (perturbation-
                                    H0  B 
                         d e                         weighted) initial conditions
                   B (0) A(t )
                                                                                                9


    Linear Response Theory: Dynamic 2.

 Now consider a more general time-dependent perturbation Fe(t)
 Simplest general form of linear response
                    t
                                                       Value at time t is a sum of the
         A(t )        dt  AB (t  t ) Fe (t )   responses to the perturbation over the
                    
                                                       entire history of the system

 For the protocol previously discussed (shut off field at t = 0)
                            0
            A(t )         dt AB (t  t)
                           
                           
                       d AB ( )
                           t
   • thus     

               d AB ( )          B(0) A(t )          AB (t )    B(0) A(t )
              t
                                                                                                    10
                                                                      t
Perturbation-Response Protocols                            A(t )     dt        B(0) A(t ) Fe (t )
                                                                      
 Turn on perturbation at t = 0, and keep
  constant thereafter
   • measured response is proportional to integral
     of time-integrated correlation function                               0
                                                                           t 
                                                         A(t )                 dt  B(0) A(t )
 Apply as -function pulse at t = 0,                                          0

  subsequent evolution proceeding normally
   • measured response proportional to time
     correlation function itself                                           0

                                                            A(t )    B(0) A(t )
 Use a sinusoidally oscillating perturbation
   • measured response proportional to Fourier-
     Laplace transformed correlation functions at
     the applied frequency                                        0
                                                              t 
   • extrapolate results from several frequencies
                                                  A(t )     dt eit B(0) A(t )
     to zero-frequency limit                                    0
                                                                                         11


                         Synthetic NEMD
 Perturb usual equations of motion in some way
    • Artificial “synthetic” perturbation need not exist in nature
 For transport coefficient of interest Lij, Ji = LijXj
    • Identify the Green-Kubo relation for the transport coefficient
                                                                
               Lij   J i ( ) J j (0) d             e.g., D   vx ( )  vx (0) d
                     0                                           0
    • Invent a fictitious field Fe, and its coupling to the system such that
      the dissipative flux is Jj
                                     H 0   J j Fe
                                       ad

    • ensure that
        equations of motion correspond to an incompressible phase space
        equations of motion are consistent with periodic boundaries
        equations of motion do not introduce inhomogeneities
    • apply a thermostat
    • couple Fe to the system and compute the steady-state average J i (t )
    • then                  J i (t )
              Lij  lim lim
                   Fe  0 t     Fe
                                                                        12



                            Phase Space
 Underlying development assumes that equations of motion
  correspond to an incompressible phase space
              q  q  p  p  0
 This can be ensured by having the perturbation derivable
  from a Hamiltonian
          H ne  H  A (p, q)  f (t )

             H ne
          q        p / m  Ap  f (t )
              p
              H ne
          p        F(q)  Aq  f (t )
               q
 Most often the equations of motion are not derivable from a
  Hamiltonian
   • but are still formulated to be compatible with an incompressible
     phase space
                                                                 13




  Diffusion: An Inhomogeneous Approach

  Artificially distinguish particles by “color”
  Introduce a species-changing plane




Molecules moving this way
across wall get colored red   Those crossing this way get blue
                                                                14




Diffusion: An Inhomogeneous Approach

 Artificially distinguish particles by “color”
 Introduce a species-changing plane




                        Considering periodic boundaries,
                        this creates a color gradient
 Problems
   • Difficult to know form of inhomogeneity in color profile
   • Cannot be extended to multicomponent diffusion
                                                                       15


              Self-Diffusion: Perturbation
 Green-Kubo relation
                                    
     D   vx ( )  vx (0) d   rx ( )  vx (0) d
          0                           0

 Label each molecule with one of two “colors”
   • each color given to half the molecules
                                                               f
 Apply Hamiltonian perturbation                                   f
                      N
        H  H 0   ci rix f (t )
                      i 1

 New equations of motion
        q  p/m
                                pix  Fix  ci f (t )
                               
        p  F(q)  Aq  f (t ) 
                                pi  y , z   Fi  y , z 
                               

 System remains homogeneous
                                                                                           16


                Self-Diffusion: Response
 Appropriate response variable is the “color current”
                         1 N
               J x (t )   ci vix (t )
                         V i 1

 According to linear response theory
                       t
       J x (t )  V  ds J x (t  s) J x (0)   0
                                                    f (s)                       f
                       0                                                               f
 In the canonical ensemble
                            1
                             2  i j xi
      J x (t ) J x (0)           c c v (t )vxj (0)
                           V i, j
                            1
                              ci2 vxi (t )vxi (0)
                           V2 i
                           N
                             vx (t )vx (0)
                           V2
                                                             1              J x (t )
                                                        D       lim lim
 Back to Green-Kubo relation                                 t  f 0     F
                                                                   17

                           Thermostatting
 External field does work on the system
   • this must be dissipated to reach steady state
 Thermostat based on velocity relative to total current density
   • “peculiar velocity”
         pix  pix  ci
         ˆ                 1
                          Nm    c j p jx
             pix  ci J x / m 

   • constrain kinetic energy
        p2 / m  3NkT
         ˆ

   • modified equations of motion
        qi  pi / m
        pi  Fi  e xci f   pi
                              ˆ

   • thermostatting multiplier

         mFi  pi
                   ˆ
           pi  pi
                 ˆ
                      Shear Viscosity:                               18




                 Boundary-Driven Algorithm
 Homogeneous algorithm for boundary-driven shear is possible
   • unique to shear viscosity
 Lees-Edwards shearing periodic boundaries (sliding brick)
   • Image cells in plane above and below central cell move

                                             dv
   • Image velocity given by shear rate   dy
                                              x




   • Peculiar velocity of all images equal                      Ly   L
           pix  pix   Ly
           ˆ




                                         Ly   L
                      Shear Viscosity:                          19




                 Boundary-Driven Algorithm
 Homogeneous algorithm for boundary-driven shear is possible
   • unique to shear viscosity
 Lees-Edwards shearing periodic boundaries (sliding brick)
   • Image cells in plane above and below central cell move

                                             dv
   • Image velocity given by shear rate   dy
                                              x




   • Peculiar velocity of all images equal
           pix  pix   Ly
           ˆ
                      Shear Viscosity:                          20




                 Boundary-Driven Algorithm
 Homogeneous algorithm for boundary-driven shear is possible
   • unique to shear viscosity
 Lees-Edwards shearing periodic boundaries (sliding brick)
   • Image cells in plane above and below central cell move

                                             dv
   • Image velocity given by shear rate   dy
                                              x




   • Peculiar velocity of all images equal
           pix  pix   Ly
           ˆ
                      Shear Viscosity:                          21




                 Boundary-Driven Algorithm
 Homogeneous algorithm for boundary-driven shear is possible
   • unique to shear viscosity
 Lees-Edwards shearing periodic boundaries (sliding brick)
   • Image cells in plane above and below central cell move

                                             dv
   • Image velocity given by shear rate   dy
                                              x




   • Peculiar velocity of all images equal
           pix  pix   Ly
           ˆ

 Try the applet
                                                22




   Lees-Edwards Boundary Conditions


                                          dvx
                           L        
                                          dy




                          Molecule exiting
                          here, in middle of
                          central cell
 L
                                                 23




   Lees-Edwards Boundary Conditions


                                           dvx
                            L        
                                           dy




                          Is replaced by one
                          here, shifted over
                          toward the edge of
 L                      the cell
                       Shift distance = Lt
                                                24




   Lees-Edwards Boundary Conditions


                                          dvx
                           L        
                                          dy




                          And with a velocity
                          that is modified
                          according to the
 L                      shear rate
                                                                                                        25




         Lees-Edwards Boundary: API

                           User's Perspective on the Molecular Simulation API


                                                Simulation


Space   Controller                      Phase            Species         Potential   Display   Device


        Integrator   MeterAbstract    Boundary        Configuration
                                                                                              26




           Lees-Edwards Boundary: Java Code
public class Space2D.BoundarySlidingBrick extends
       Space2D.BoundaryPeriodicSquare
public void nearestImage(Vector dr) {
  double delrx = delvx*timer.currentValue();
  double cory;
  cory = (dr.y > 0.0) ? Math.floor(dr.y/dimensions.y+0.5):Math.ceil(dr.y/dimensions.y-0.5);
  dr.x -= cory*delrx;
  dr.x -= dimensions.x * ((dr.x > 0.0) ? Math.floor(dr.x/dimensions.x+0.5) :
                           Math.ceil(dr.x/dimensions.x-0.5));
  dr.y -= dimensions.y * cory;
}


public void centralImage(Coordinate c) {
  Vector r = c.r;
  double cory = (r.y > 0.0) ? Math.floor(r.y/dimensions.y) : Math.ceil(r.y/dimensions.y-1.0);
  double corx = (r.x > 0.0) ? Math.floor(r.x/dimensions.x) : Math.ceil(r.x/dimensions.x-1.0);
  if(corx==0.0 && cory==0.0) return;
  double delrx = delvx*timer.currentValue();
  Vector p = c.p;
  r.x -= cory*delrx;
  r.x -= dimensions.x * corx;
  r.y -= dimensions.y * cory;
  p.x -= cory*delvx;
}
                                                          27




  Limitations of Boundary-Driven Shear

 No external field in equations of motion
   • cannot employ response theory to link to viscosity
 Lag time in response of system to initiation of shear
   • cannot be used to examine time-dependent flows
 A fictitious-force method is preferable
                                                               28


DOLLS-Tensor Hamiltonian: Perturbation
 An arbitrary fictitious shear field can be imposed via the
  DOLLS-tensor Hamiltonian
                     N
         H  H 0   qipi   u(t ) 
                                      T
                          
                    i 1

 Equations of motion
         qi  pi / m  qi  u
         pi  Fi  u  pi

   • must be implemented with compatible PBC
 Example: Simple Couette shear

            0 0 0              qi  pi / m   qiy e x
      u    0 0              pi  Fi   pixe y
                  
            0 0 0
                  
                                                         29




 DOLLS-Tensor Hamiltonian: Response
 Appropriate response variable is the pressure tensor
                1 N 1          N
        P (t )   m pipi  2  rij Fij
                            1
                V i 1        i, j

 According to linear response theory
                   t
      P(t )   V  ds P(t  s)P(0) 0  u( s)
                                       
                   0

 Shear viscosity, via Green-Kubo
                              Pxy (t )
             lim lim
                t    0      
                                                                                     30


                  SLLOD Formulation
 DOLLS-tensor formulation fails in more complex situations
   • non-linear regime
   • evaluation of normal-stress differences
   • a simple change fixes things up
 SLLOD Equations of motion                                    DOLLS
         qi  pi / m  qi  u                             qi  pi / m  qi  u
         pi  Fi  pi  u                                 pi  Fi  u  pi
                             Only change
 Example: Simple Couette shear
            0 0 0              qi  pi / m   qiy e x   qi  pi / m   qiy e x
      u    0 0              pi  Fi   piy e x
                                                         pi  Fi   pixe y
            0 0 0
                  
 Methods equivalent for irrotational flows                u   u 
                                                                        T
                                                                      31


                         Application

 NEMD usually introduces exceptionally large strain rates
   • 108 sec-1 or greater
                                          
                                                1/ 2
                                         ms 2
   • dimensionless strain rate  *      
   • thus, e.g.,
       m = 30g/mol; s = 3A; /k = 100K; * = 1.0   = 51011 sec-1
 Shear-thinning observed even in simple fluids at these rates
 Very important to extrapolate to zero shear
             


Newtonian

                                                        *1/ 2

				
DOCUMENT INFO
Shared By:
Categories:
Tags:
Stats:
views:6
posted:3/23/2011
language:French
pages:31