## 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