UCL WIKI

UCL Logo
Child pages
  • Exercise - Molecular dynamics
Skip to end of metadata
Go to start of metadata

Here is an MD program for the dynamics of two planets 

The Verlet MD code for planets
!****************************************************************************
!
!  PROGRAM: MD_verlet
!
!****************************************************************************
    program MD_verlet
    implicit none
    ! Variables
    !
    integer,parameter :: N=2,dim=3
    !
    real(8),parameter :: au = 4.503e9 ! m 
    real(8),parameter :: msun = 1.9891e30
    real(8),parameter :: mearth = 3.0e-6*msun  ! suns
    real(8),parameter :: G = 6.67e-11 ! N m^2/kg^2
    real(8),parameter :: dt = 10.0
    integer,parameter :: Ntime = 10000
    !
    real(8) :: r(N,dim),v(N,dim),a(N,dim),force(dim),dist,m(N)
    integer :: i,j,itime
    ! Body of MD_verlet
    ! intialization
    !
    r(1,:) = 0 
    v(1,:) = 0
    a(1,:) = 0
    m(1) = msun
    !
    r(2,:) = (/0.0d0,0.4d0,0.0d0/)*au
    v(2,:) = (/200250.0d0,0.0d0,0.0d0/)
    a(2,:) = 0
    m(2) = 0.055*MEarth
    !
    !U = -G*M*m/r
    !
    do itime = 1,Ntime
      !
      if ( mod(itime,10)==0 ) write(8,"(i4)") N
      if ( mod(itime,10)==0 ) write(8,*)
      !
      do i = 1,N
          !
          force = 0
          do j = 1,N
            !
            if (i==j) cycle
            !
            dist = sqrt( sum(  ( r(i,:)-r(j,:) )**2 ) )
            force(:) = force(:) - G*m(i)*m(j)/dist**2*(r(i,:)-r(j,:))/dist
            !
          enddo
          !
          r(i,:) = r(i,:) + v(i,:)*dt + 0.5d0*dt*dt*a(i,:)
          v(i,:) = v(i,:) + 0.5d0*dt*(force(:)/m(i) + a(i,:))
          a(i,:) = force(:)/m(i)
          !
          if ( mod(itime,10)==0 ) write(8,"('C ',<dim>(1x,e14.7))") r(i,:)/au*10.0
          !
      enddo
      !
    enddo
    !
    end program MD_verlet

 

Here is an example (MD_verlet_atoms.f90):

MD_verlet_atoms.f90
!****************************************************************************
!
!  PROGRAM: MD_verlet
!
!  PURPOSE:  Entry point for the console application.
!
!****************************************************************************
    program MD_verlet_atoms
    implicit none
    ! Variables
    !
    integer,parameter :: N=4,dim=3
    !
    real(8),parameter :: au = 4.503e9 ! m 
    real(8),parameter :: msun = 1.9891e30
    real(8),parameter :: eps = 0.1
    real(8),parameter :: sigma  = 1.4
    real(8),parameter :: Na = 1.0 ! 6.0221415e+23 ! 
    real(8),parameter :: dt = 0.01
    integer,parameter :: Ntime = 100000
    !
    real(8) :: r(N,dim),v(N,dim),a(N,dim),force(dim),dist,m(N)
    integer :: i,j,itime
    ! Body of MD_verlet
    ! intialization
    !
    r(1,:) = 0 
    v(1,:) = 0
    a(1,:) = 0
    m(1) = 12.0/Na
    !
    r(2,:) = (/0.0d0,1.2d0,0.5d0/)
    v(2,:) = 0
    a(2,:) = 0
    m(2) = 12.0/Na
    !
    r(3,:) = (/0.0d0,2.4d0,-0.5d0/)
    v(3,:) = 0
    a(3,:) = 0
    m(3) = 12.0/Na
    !
    r(4,:) = (/0.0d0,3.6d0,0.5d0/)
    v(4,:) = 0
    a(4,:) = 0
    m(4) = 12.0/Na
    !
    do itime = 1,Ntime
      !
      if ( mod(itime,10)==0 ) write(8,"(i4)") N
      if ( mod(itime,10)==0 ) write(8,*)
      !
      do i = 1,N
          !
          force = 0
          do j = 1,N
            !
            if (i==j) cycle
            !
            dist = sqrt( sum(  ( r(i,:)-r(j,:) )**2 ) )
            force(:) = force(:) + 0.5d0*eps*(-12.0*sigma**12/dist**13+2.0*6.0*sigma**6/dist**7)*(r(j,:)-r(i,:))/dist
            !
          enddo
          !
          r(i,:) = r(i,:) + v(i,:)*dt + 0.5d0*dt*dt*a(i,:)
          v(i,:) = v(i,:) + 0.5d0*dt*(force(:)/m(i) + a(i,:))
          a(i,:) = force(:)/m(i)
          !
          if ( mod(itime,10)==0 ) write(8,"('C ',<dim>(1x,e14.7))") r(i,:)
          !
      enddo
      !
    enddo
    !
    end program MD_verlet_atoms


Exercise 1: modify the MD programs to include more planets or atoms (see program below)


Redirect the output to the file MD.out (see program atributed)

> MD.out

and plot the potential energy values using a gnuplot script MD.gpl (see https://wiki.ucl.ac.uk/display/GradPrgCrs/gnuplot+scripts)

gnuplot –persist MD.gpl

 

Apart from the output file MD.out you should obtain the file “fort.8” with Cartesian coordinates of 500 particles listed as follows:

fort.8
   2
 
C   0.0000000E+00  0.0000000E+00  0.0000000E+00
C   0.2019562-307  0.4000000E+01  0.0000000E+00
   2
 
C   0.0000000E+00  0.0000000E+00  0.0000000E+00
C   0.1998854-306  0.4000000E+01  0.0000000E+00
   2
 
C   0.0000000E+00  0.0000000E+00  0.0000000E+00
C   0.3795752-306  0.4000000E+01  0.0000000E+00
   2
 
 
Use wxmacmolplt o visualize the molecular dynamics:
wxmacmolplt   fort.8

or with molden

/nfs/workspaces/ucapsyu/Fortran95/sergey/Verlet/Verlet/molden fort.8

 

An object-oriented version of this code consists of three parts MD_verlet.f90driver.f90constants.f90.

Exercise 3: Open these files with Code::Blocks, compile and run. 

MD_verlet.f90
!****************************************************************************
!
!  PROGRAM: MD_verlet
!
!****************************************************************************
 module MD_verlet
   !
   use constants
   !
   public Verlet
   !
   integer,parameter :: dim=3   ! dimension of the problem
   real(8),parameter :: dt = 10.0 ! time step 
   integer,parameter :: Ntime = 10000  ! Number of time steps
   !
   ! declare a new typo to describe atoms/planets
   type objectT
     !
     real(8) :: r(3) ! position 
     real(8) :: v(3) ! velocity 
     real(8) :: a(3) ! acceleration
     real(8) :: mass ! 
   contains
     procedure :: Evolution
   end type objectT
   !
   contains 
   !
   subroutine Verlet 
     !
     implicit none
     !
     ! Variables
     !
     integer,parameter :: N=2
     !
     real(8) :: force(dim),dist
     integer :: i,j,itime
     !
     type(objectT) :: object(2)
   
     ! Body of MD_verlet
   
     ! intialization
     !
     call constantsInitialize
     !
     object(1)%r(:) = 0 
     object(1)%v(:) = 0 
     object(1)%a(:) = 0 
     object(1)%mass = msun
     !
     object(2)%r = (/0.0d0,0.4d0,0.0d0/)*au
     object(2)%v = (/200250.0d0,0.0d0,0.0d0/)
     object(2)%a = 0
     object(2)%mass = 0.055*MEarth
     !
     !U = -G*M*m/r
     !
     do itime = 1,Ntime
       !
       if ( mod(itime,10)==0 ) write(8,"(i4)") N
       if ( mod(itime,10)==0 ) write(8,*)
       !
       do i = 1,N
           !
           force = 0
           do j = 1,N
             !
             if (i==j) cycle
             !
             dist = sqrt( sum(  ( object(i)%r(:)-object(j)%r(:) )**2 ) )
             force(:) = force(:) - G*object(i)%mass*object(j)%mass/dist**2*(object(i)%r(:)-object(j)%r(:))/dist
             !
           enddo
           !
           call object(i)%Evolution(force)
           !
           if ( mod(itime,10)==0 ) write(8,"('C ',<dim>(1x,e14.7))") object(i)%r(:)/au*10.0
           !
       enddo
       !
     enddo
     !
   end subroutine Verlet
   !
   subroutine Evolution(object,force)
      !
      implicit none
      !
      class(objectT), intent(inout) :: object
      real(8) :: force(dim),dt
      integer :: getAge
      !
      object%r(:) = object%r(:) + object%v(:)*dt + 0.5d0*dt*dt*object%a(:)
      object%v(:) = object%v(:) + 0.5d0*dt*(force(:)/object%mass + object%a(:) )
      object%a(:) = force(:)/object%mass
      !
   end subroutine
   !
 end module MD_verlet
constants.f90
module constants
  implicit none
  private
  public constantsInitialize
  public planck,avogno,vellgt,boltz,bohr,todebye
  public hartree,ev, msun, au, mearth, G
  !
  real(8),parameter :: au = 4.503e9 ! m 
  real(8),parameter :: msun = 1.9891e30
  real(8),parameter :: mearth = 3.0e-6*msun  ! suns
  real(8),parameter :: G = 6.67e-11 ! N m^2/kg^2
  ! physical constants -- All constants updated 21 March 2012 from the NIST
  ! website http://physics.nist.gov/cuu/Constants/index.html
  real(8), parameter :: planck     =  6.62606957e-27           ! Planck constant in (non-SI) erg*second
  real(8), parameter :: avogno     =  6.02214129e+23           ! Avogadro constant
  real(8), parameter :: vellgt     =  2.99792458e+10           ! Speed of light constant in (non-SI) cm/second
  real(8), parameter :: boltz      =  1.3806488e-16            ! Boltzmann constant in (non-SI) erg/Kelvin
  real(8), parameter :: bohr       =  0.52917720859d0          ! bohr constant in Angstroms
  real(8), parameter :: hartree    =  219474.6313705d0         ! hartree in cm-1
  real(8), parameter :: ev         =  8065.54465d0             ! ev in cm-1
  real(8), parameter :: uma        =  1.660538921e-24          ! unified atomic mass unit [=mass(C-12)/12 ] in grams
  real(8), parameter :: todebye    =  2.541765d0               ! a.u. in debye
  
  contains
    subroutine constantsInitialize  ! This initialization routine is now superfluous but kept for compatibility
     !
    end subroutine constantsInitialize
end module constants
driver.f90
 program driver
 use MD_verlet
  !
  ! Molecular dynamics, Verlet
  call Verlet
  !
 end program driver 


 

Here is an example of an MD simulation visualized with molden: 

 

 

 

  • No labels