Molecular Dynamics Algorithms

Revision as of 09:00, 23 July 2021 by Kulhanek (talk | contribs)
(diff) ← Older revision | Latest revision (diff) | Newer revision → (diff)

Navigation: Documentation / Technical Details / Molecular Dynamics Algorithms


AMBER

Both sander and pmemd from the AMBER package use the leap-frog integration algorithm. The algorithm is demonstrated for the SANDER driver:

  1. Initialization (sander.F90)
  2. call pmf_sander_init_preinit(mdin,natom,nres,ntb,nstlim,dt,temp0,a,b,c,alpha,beta,gamma)
    ! topology population
    do i=1,nres
        call pmf_sander_set_residue(i,ih(m02-1+i),ix(i02-1+i))
    end do
    do i=1,natom
        call pmf_sander_set_atom(i,ih(m04-1+i),ih(m06-1+i))
    end do
    call pmf_sander_finalize_preinit(natom,x(lmass),x(lcrd))
    call pmf_sander_cst_init_collisions(ntc,nbonh,ix(iifstwt),ix(iibh),ix(ijbh),x(l50))
    call pmf_sander_init(natom,x(lmass),x(lcrd))
    
  3. MD Loop (runmd.F90)
    1. Calculate Forces
    2. call force(x,f,ener)  ! AMBER MM potential and forces
      ! x - in t
      ! v - in t-dt/2
      ! f - in t
      ! ener%pot%tot - in t
      ! ener%kin%tot - in t-dt
      if( ntb .ne. 0 ) then
          call pmf_sander_update_box(a,b,c,alpha,beta,gamma)
      end if
      call pmf_sander_force(natom,x,v,f,ener%pot%tot,ener%kin%tot,pmfene)
      ener%pot%tot = ener%pot%tot + pmfene
      
    3. Update Velocities + thermostat
    4. ! backup old velocities
      vold(t-dt/2) = v(t-dt/2)
      ! update velocities
      v(t+dt/2) = v(t-dt/2) + a(t)*dt
      
    5. Update Positions + barostat
    6. ! backup positions
      f(t) = x(t)
      ! update positions
      x(t+dt) = x(t) + v(t+dt/2)*dt
      
    7. Apply Constraints for SHAKE or PMFLib constraints
    8. call pmf_sander_constraints(natom,x,con_modified)
      call shake
      ! f - old positions in t
      ! x - updated positions in t+dt
      ! re-evaluate velocities
      v(t+dt/2) = [x(t+dt) - f(t)]/dt
      
    9. Calculate Kinetic Energy
    10. v(t) = [v(t+dt/2) + vold(t-dt/2)] / 2
      ener%kin%tot <- v(t)
      
  4. Finalization (sander.F90)
  5. call pmf_sander_finalize()