# UCL WIKI

##### Child pages
• Exercise - Molecular dynamics
Go to start of banner

# Exercise - Molecular dynamics

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