MD Program

This is a simple Molecular Dynamics program.

Program

program main

*****************************************************************************80

! MAIN is the main program for MD.

Discussion:

MD implements a simple molecular dynamics simulation.

The velocity Verlet time integration scheme is used.

The particles interact with a central pair potential.

Based on a FORTRAN90 program by Bill Magro.

Calculating distances, forces, potential- and kinetic- energies have been broken out into subroutines to make this a more interesting example to analyze with a profiler.

2018-05-25 Oliver Stueker

Usage:

md nd np step_num dt

where

  • nd is the spatial dimension (2 or 3);
  • np is the number of particles (500, for instance);
  • step_num is the number of time steps (500, for instance).
  • dt is the time step (0.1 for instance )

Licensing:

This code is distributed under the GNU LGPL license.

Modified:

25 May 2018 by Oliver Stueker

Author:

John Burkardt
Call to:timestamp(), s_to_i4(), s_to_r8(), initialize(), update(), compute(), calc_distance(), calc_pot(), calc_force(), calc_kin(), r8mat_uniform_ab()

Subroutines and functions

subroutine compute(np, nd, pos, vel, mass, f, pot, kin)

*****************************************************************************80

! COMPUTE computes the forces and energies.

Discussion:

The computation of forces and energies is fully parallel.

The potential function V(X) is a harmonic well which smoothly saturates to a maximum value at PI/2:

v(x) = ( sin ( min ( x, PI/2 ) ) )^2

The derivative of the potential is:

dv(x) = 2.0D+00 * sin ( min ( x, PI/2 ) ) * cos ( min ( x, PI/2 ) )
= sin ( 2.0 * min ( x, PI/2 ) )

Licensing:

This code is distributed under the GNU LGPL license.

Modified:

15 July 2008

Author:

John Burkardt

Parameters:

Input, integer ( kind = 4 ) NP, the number of particles.

Input, integer ( kind = 4 ) ND, the number of spatial dimensions.

Input, real ( kind = 8 ) POS(ND,NP), the positions.

Input, real ( kind = 8 ) VEL(ND,NP), the velocities.

Input, real ( kind = 8 ) MASS, the mass.

Output, real ( kind = 8 ) F(ND,NP), the forces.

Output, real ( kind = 8 ) POT, the total potential energy.

Output, real ( kind = 8 ) KIN, the total kinetic energy.

Options:
  • np [integer,optional/default=shape(pos,1)]
  • nd [integer,optional/default=shape(pos,0)]
Parameters:
  • pos (nd,np) [real]
  • vel (nd,np) [real]
  • mass [real]
  • f (nd,np) [real]
  • pot [real]
  • kin [real]
Called from:

main
Call to:

calc_distance(), calc_pot(), calc_force(), calc_kin()
subroutine calc_pot(d2, pot)

Calculate potential energy for truncated distance D2.

Parameters:
Input, real ( kind = 8 ) D2, trucated distance. In/Out, real ( kind = 8 ) POT, potential energy.
Parameters:
  • d2 [real]
  • pot [real]
Called from:

compute(), main
subroutine calc_force(np, nd, i, d, d2, rij, f)

Add particle J’s contribution to the force on particle I.

Parameters:
Input, integer ( kind = 4 ) NP, the number of particles. Input, integer ( kind = 4 ) ND, the number of spatial dimensions. Input, integer ( kind = 4 ) I, index of particle I. Input, real ( kind = 8 ) D, distance. Input, real ( kind = 8 ) D2, trucated distance. Input, real ( kind = 8 ) Rij(nd), distance vector In/Out, real ( kind = 8 ) F(ND,NP), the forces.
Options:
  • np [integer,optional/default=shape(f,1)]
  • nd [integer,optional/default=len(rij)]
Parameters:
  • i [integer]
  • d [real]
  • d2 [real]
  • rij (nd) [real]
  • f (nd,np) [real]
Called from:

compute(), main
subroutine calc_kin(np, nd, vel, mass, kin)

Compute the total kinetic energy.

Parameters:
Input, integer ( kind = 4 ) NP, the number of particles. Input, integer ( kind = 4 ) ND, the number of spatial dimensions. Input, real ( kind = 8 ) VEL(ND,NP), the velocities. Input, real ( kind = 8 ) MASS, the mass. Output, real ( kind = 8 ) KIN, the total kinetic energy.
Options:
  • np [integer,optional/default=shape(vel,1)]
  • nd [integer,optional/default=shape(vel,0)]
Parameters:
  • vel (nd,np) [real]
  • mass [real]
  • kin [real]
Called from:

compute(), main
subroutine initialize(np, nd, pos, vel, acc)
Options:
  • np [integer,optional/default=shape(pos,1)]
  • nd [integer,optional/default=shape(pos,0)]
Parameters:
  • pos (nd,np) [real]
  • vel (nd,np) [real]
  • acc (nd,np) [real]
Called from:

main
Call to:

r8mat_uniform_ab()
subroutine r8mat_uniform_ab(m, n, a, b, seed, r)
Options:
  • m [integer,optional/default=shape(r,0)]
  • n [integer,optional/default=shape(r,1)]
Parameters:
  • a [real]
  • b [real]
  • seed [integer]
  • r (m,n) [real] ::

    ’
Called from:

initialize(), main
subroutine s_to_i4(s, value, ierror, length)
Parameters:
  • s [character] ::

    ’

  • value [integer]
  • ierror [integer]
  • length [integer]
Called from:

main
subroutine s_to_r8(s, r8)
Parameters:
  • s [character] ::

    ’

  • r8 [real]
Called from:

main
subroutine timestamp()
Called from:main
subroutine update(np, nd, pos, vel, f, acc, mass, dt)
Options:
  • np [integer,optional/default=shape(pos,1)]
  • nd [integer,optional/default=shape(pos,0)]
Parameters:
  • pos (nd,np) [real]
  • vel (nd,np) [real]
  • f (nd,np) [real]
  • acc (nd,np) [real]
  • mass [real]
  • dt [real]
Called from:

main