## MD is the main program for the molecular dynamics simulation.
#
#  Discussion:
#    MD implements a simple molecular dynamics simulation.
#
#    The velocity Verlet time integration scheme is used.
#    The particles interact with a central pair potential.
#
#  Licensing:
#    This code is distributed under the GNU LGPL license.
#  Modified:
#    26 December 2014
#
#  Author:
#    John Burkardt
#
#  Parameters:
#    Input, integer D_NUM, the spatial dimension.
#    A value of 2 or 3 is usual.
#    The default value is 3.
#
#    Input, integer P_NUM, the number of particles.
#    A value of 1000 or 2000 is small but "reasonable".
#    The default value is 500.
#
#    Input, integer STEP_NUM, the number of time steps.
#    A value of 500 is a small but reasonable value.
#    The default value is 500.
#
#    Input, real DT, the time step.
#    A value of 0.1 is large; the system will begin to move quickly but the
#    results will be less accurate.
#    A value of 0.0001 is small, but the results will be more accurate.
#    The default value is 0.1.
#
import platform
from time import clock
import numpy as np
from sys import exit
import time
def timestamp ( ):
    t = time.time ( )
    print ( time.ctime ( t ) )
    return None
# end def
# ===================
# Main section of code
# ====================
timestamp()
print('')
print('MD_TEST')
print('  Python version: %s' % (platform.python_version( ) ))
print('  Test the MD molecular dynamics program.')
# Initialize variables
d_num = 3
p_num = 500
step_num = 50
dt = 0.1
mass = 1.0
rmass = 1.0 / mass
wtime1 = clock( )
#  output:
print('' )
print('MD' )
print('  Python version: %s' % (platform.python_version( ) ) )
print('  A molecular dynamics program.' )
print('' )
print('  D_NUM, the spatial dimension, is %d' % (d_num) )
print('  P_NUM, the number of particles in the simulation is %d.' % (p_num) )
print('  STEP_NUM, the number of time steps, is %d.' % (step_num) )
print('  DT, the time step size, is %g seconds.' % (dt) )
print('' )
print('  At each step, we report the potential and kinetic energies.' )
print('  The sum of these energies should be a constant.' )
print('  As an accuracy check, we also print the relative error' )
print('  in the total energy.' )
print('' )
print('      Step      Potential       Kinetic        (P+K-E0)/E0' )
print('                Energy P        Energy K       Relative Energy Error')
print('')
step_print_index = 0
step_print_num = 10
step_print = 0
for step in range(0, step_num+1):
    if (step == 0):
        # Initialize
        #  Positions
        seed = 123456789
        a = 0.0
        b = 10.0
        i4_huge = 2147483647
        if (seed < 0):
            seed = seed + i4_huge
        # end if
        if (seed == 0):
             print('' )
             print( 'R8MAT_UNIFORM_AB - Fatal error!')
             print('  Input SEED = 0!' )
             sys.ext('R8MAT_UNIFORM_AB - Fatal error!')
        # end if
        pos = np.zeros( (d_num, p_num) )
        for j in range(0, p_num):
            for i in range(0, d_num):
                k = (seed // 127773)
                seed = 16807 * (seed - k * 127773) - k * 2836
                seed = (seed % i4_huge)
                if (seed < 0):
                       seed = seed + i4_huge
                # end if
                pos[i,j] = a + (b - a) * seed * 4.656612875E-10
            # end for
        # end for
        #  Velocities
        vel = np.zeros([ d_num, p_num ])
        #  Accelerations
        acc = np.zeros([ d_num, p_num ])
    else:
        # Update
        #  Update positions
        for i in range(0, d_num):
             for j in range(0, p_num):
                pos[i,j] = pos[i,j] + vel[i,j] * dt + 0.5 * acc[i,j] * dt * dt
            # end for
        # end for
        #  Update velocities
        for i in range(0, d_num):
            for j in range(0, p_num):
                 vel[i,j] = vel[i,j] + 0.5 * dt * ( force[i,j] * rmass + acc[i,j] )
            # end for
        # end for
         #  Update accelerations.
        for i in range(0, d_num):
             for j in range(0, p_num):
                acc[i,j] = force[i,j] * rmass
            # end for
        # end for
    # endif
    # compute force, potential, kinetic
    force = np.zeros([ d_num, p_num ])
    rij = np.zeros(d_num)
    potential = 0.0
    for i in range(0, p_num):
        #  Compute the potential energy and forces.
        for j in range(0, p_num):
            if (i != j):
                #  Compute RIJ, the displacement vector.
                for k in range(0, d_num):
                    rij[k] = pos[k,i] - pos[k,j]
                # end for
                #  Compute D and D2, a distance and a truncated distance.
                d = 0.0
                for k in range(0, d_num):
                    d = d + rij[k] ** 2
                # end for
                d = np.sqrt(d)
                d2 = min(d, np.pi / 2.0)
                #  Attribute half of the total potential energy to particle J.
                potential = potential + 0.5 * np.sin(d2) * np.sin(d2)
                #  Add particle J's contribution to the force on particle I.
                for k in range(0, d_num):
                    force[k,i] = force[k,i] - rij[k] * np.sin(2.0 * d2) / d
                # end for
            # end if
        # end for
    # end for
    #  Compute the kinetic energy
    kinetic = 0.0
    for k in range(0, d_num):
        for j in range(0, p_num):
            kinetic = kinetic + vel[k,j] ** 2
        # end for
     # end for
    kinetic = 0.5 * mass * kinetic
    if (step == 0):
         e0 = potential + kinetic
    # endif
    if (step == step_print):
        rel = (potential + kinetic - e0) / e0
        print('  %8d  %14f  %14f  %14g' % (step, potential, kinetic, rel) )
        step_print_index = step_print_index + 1
        step_print = (step_print_index * step_num) // step_print_num
    #end if
# end step
wtime2 = clock( )
print('')
print('    Elapsed wall clock time = %g seconds.' % (wtime2 - wtime1) )
#  Terminate
print('')
print('MD_TEST')
print('  Normal end of execution.')
timestamp ( )
# end if