Academic Integrity: tutoring, explanations, and feedback — we don’t complete graded work or submit on a student’s behalf.

function md ( nd, np, step_num, dt ) %******************************************

ID: 3778932 • Letter: F

Question

function md ( nd, np, step_num, dt )

%*****************************************************************************80
%
%% 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.
%
% Based on a FORTRAN90 program by Bill Magro.
%
% 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:
%
% 27 December 2014
%
% Author:
%
% John Burkardt
%
% Parameters:
%
% Input, integer ND, the spatial dimension.
% A value of 2 or 3 is usual.
% The default value is 3.
%
% Input, integer NP, 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.
%
if ( nargin < 1 )
nd = 3;
elseif ( ischar ( nd ) )
nd = str2num ( nd );
end

if ( nargin < 2 )
np = 500;
elseif ( ischar ( np ) )
np = str2num ( np );
end

if ( nargin < 3 )
step_num = 500;
elseif ( ischar ( step_num ) )
step_num = str2num ( step_num );
end

if ( nargin < 4 )
dt = 0.1;
elseif ( ischar ( dt ) )
dt = str2num ( dt );
end

mass = 1.0;
%
% Report to the user.
%
timestamp ( );
fprintf ( 1, ' ' );
fprintf ( 1, 'MD ' );
fprintf ( 1, ' MATLAB version ' );
fprintf ( 1, ' A molecular dynamics program. ' );
fprintf ( 1, ' ' );
fprintf ( 1, ' ND, the spatial dimension, is %d ', nd );
fprintf ( 1, ' NP, the number of particles in the simulation is %d. ', np );
fprintf ( 1, ' STEP_NUM, the number of time steps, is %d. ', step_num );
fprintf ( 1, ' DT, the time step size, is %g seconds. ', dt );
%
% This is the main time stepping loop:
% Compute forces and energies,
% Update positions, velocities, accelerations.
%
fprintf ( 1, ' ' );
fprintf ( 1, ' At each step, we report the potential and kinetic energies. ' );
fprintf ( 1, ' The sum of these energies should be a constant. ' );
fprintf ( 1, ' As an accuracy check, we also print the relative error ' );
fprintf ( 1, ' in the total energy. ' );
fprintf ( 1, ' ' );
fprintf ( 1, ' Step Potential Kinetic (P+K-E0)/E0 ' );
fprintf ( 1, ' Energy P Energy K Relative Energy Error ' );
fprintf ( 1, ' ' );

step_print_index = 0;
step_print_num = 10;
step_print = 0;

tic;

for step = 0 : step_num

if ( step == 0 )
[ pos, vel, acc ] = initialize ( np, nd );
else
[ pos, vel, acc ] = update ( np, nd, pos, vel, force, acc, mass, dt );
end

[ force, potential, kinetic ] = compute ( np, nd, pos, vel, mass );

if ( step == 0 )
e0 = potential + kinetic;
end

if ( step == step_print )
rel = ( potential + kinetic - e0 ) / e0;
fprintf ( 1, ' %8d %14f %14f %14e ', step, potential, kinetic, rel );
step_print_index = step_print_index + 1;
step_print = floor ( ( step_print_index * step_num ) / step_print_num );
end

end
%
% Report timing.
%
wtime = toc;
fprintf ( 1, ' ' );
fprintf ( 1, ' Wall clock time = %f seconds. ', wtime );
%
% Terminate.
%
fprintf ( 1, ' ' );
fprintf ( 1, 'MD ' );
fprintf ( 1, ' Normal end of execution. ' );
fprintf ( 1, ' ' );
timestamp ( );

return
end
function [ f, pot, kin ] = compute ( np, nd, pos, vel, mass )

%*****************************************************************************80
%
%% COMPUTE computes the forces and energies.
%
% Discussion:
%
% 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.0 * 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:
%
% 30 July 2009
%
% Author:
%
% John Burkardt and Gene Cliff.
%
% Parameters:
%
% Input, integer NP, the number of particles.
%
% Input, integer ND, the number of spatial dimensions.
%
% Input, real POS(ND,NP), the positions.
%
% Input, real VEL(ND,NP), the velocities.
%
% Input, real MASS, the mass of each particle.
%
% Output, real F(ND,NP), the forces.
%
% Output, real POT, the total potential energy.
%
% Output, real KIN, the total kinetic energy.
%
f = zeros ( nd, np );

pot = 0.0;

for i = 1 : np
%
% Compute the potential energy and forces.
%
for j = 1 : np

if ( i ~= j )

for k = 1 : nd
rij(k) = pos(k,i) - pos(k,j);
end

d = 0.0;
for k = 1 : nd
d = d + rij(k)^2;
end
d = sqrt ( d );
%
% Truncate the distance.
%
d2 = min ( d, pi / 2.0 );
%
% Attribute half of the total potential energy to particle J.
%
pot = pot + 0.5 * sin ( d2 ) * sin ( d2 );
%
% Add particle J's contribution to the force on particle I.
%
for k = 1 : nd
f(k,i) = f(k,i) - rij(k) * sin ( 2.0 * d2 ) / d;
end

end

end

end
%
% Compute the total kinetic energy.
%
kin = 0.5 * mass * sum ( sum ( vel(1:nd,1:np).^2 ) );

return
end
function [ pos, vel, acc ] = initialize ( np, nd )

%*****************************************************************************80
%
%% INITIALIZE initializes the positions, velocities, and accelerations.
%
% Licensing:
%
% This code is distributed under the GNU LGPL license.
%
% Modified:
%
% 26 December 2014
%
% Author:
%
% John Burkardt.
%
% Parameters:
%
% Input, integer NP, the number of particles.
%
% Input, integer ND, the number of spatial dimensions.
%
% Output, real POS(ND,NP), the positions.
%
% Output, real VEL(ND,NP), the velocities.
%
% Output, real ACC(ND,NP), the accelerations.
%
seed = 123456789;
%
% Positions.
%
[ pos, seed ] = r8mat_uniform_ab ( nd, np, 0.0, 10.0, seed );
%
% Velocities.
%
vel = zeros ( nd, np );
%
% Accelerations.
%
acc = zeros ( nd, np );

return
end
function [ r, seed ] = r8mat_uniform_ab ( m, n, a, b, seed )

%*****************************************************************************80
%
%% R8MAT_UNIFORM_AB returns a scaled pseudorandom R8MAT.
%
% Discussion:
%
% An R8MAT is an array of R8's.
%
% Licensing:
%
% This code is distributed under the GNU LGPL license.
%
% Modified:
%
% 21 September 2006
%
% Author:
%
% John Burkardt
%
% Reference:
%
% Paul Bratley, Bennett Fox, Linus Schrage,
% A Guide to Simulation,
% Second Edition,
% Springer, 1987,
% ISBN: 0387964673,
% LC: QA76.9.C65.B73.
%
% Bennett Fox,
% Algorithm 647:
% Implementation and Relative Efficiency of Quasirandom
% Sequence Generators,
% ACM Transactions on Mathematical Software,
% Volume 12, Number 4, December 1986, pages 362-376.
%
% Pierre L'Ecuyer,
% Random Number Generation,
% in Handbook of Simulation,
% edited by Jerry Banks,
% Wiley, 1998,
% ISBN: 0471134031,
% LC: T57.62.H37.
%
% Peter Lewis, Allen Goodman, James Miller,
% A Pseudo-Random Number Generator for the System/360,
% IBM Systems Journal,
% Volume 8, Number 2, 1969, pages 136-143.
%
% Parameters:
%
% Input, integer M, N, the number of rows and columns in the array.
%
% Input, real A, B, the range of the pseudorandom values.
%
% Input, integer SEED, the integer "seed" used to generate
% the output random number.
%
% Output, real R(M,N), an array of random values between 0 and 1.
%
% Output, integer SEED, the updated seed. This would
% normally be used as the input seed on the next call.
%
r = zeros ( m, n );

i4_huge = 2147483647;

if ( seed == 0 )
fprintf ( 1, ' ' );
fprintf ( 1, 'R8MAT_UNIFORM_AB - Fatal error! ' );
fprintf ( 1, ' Input SEED = 0! ' );
error ( 'R8MAT_UNIFORM_AB - Fatal error!' );
end

for j = 1 : n
for i = 1 : m

seed = floor ( seed );

seed = mod ( seed, i4_huge );

if ( seed < 0 )
seed = seed + i4_huge;
end

k = floor ( seed / 127773 );

seed = 16807 * ( seed - k * 127773 ) - k * 2836;

if ( seed < 0 )
seed = seed + i4_huge;
end

r(i,j) = a + ( b - a ) * seed * 4.656612875E-10;

end
end

return
end
function timestamp ( )

%*****************************************************************************80
%
%% TIMESTAMP prints the current YMDHMS date as a timestamp.
%
% Licensing:
%
% This code is distributed under the GNU LGPL license.
%
% Modified:
%
% 14 February 2003
%
% Author:
%
% John Burkardt
%
t = now;
c = datevec ( t );
s = datestr ( c, 0 );
fprintf ( 1, '%s ', s );

return
end
function [ pos, vel, acc ] = update ( np, nd, pos, vel, f, acc, mass, dt )

%*****************************************************************************80
%
%% UPDATE updates positions, velocities and accelerations.
%
% Discussion:
%
% The time integration is fully parallel.
%
% A velocity Verlet algorithm is used for the updating.
%
% x(t+dt) = x(t) + v(t) * dt + 0.5 * a(t) * dt * dt
% v(t+dt) = v(t) + 0.5 * ( a(t) + a(t+dt) ) * dt
% a(t+dt) = f(t) / m
%
% Licensing:
%
% This code is distributed under the GNU LGPL license.
%
% Modified:
%
% 15 July 2008
%
% Author:
%
% John Burkardt
%
% Parameters:
%
% Input, integer NP, the number of particles.
%
% Input, integer ND, the number of spatial dimensions.
%
% Input, real POS(ND,NP), the positions.
%
% Input, real VEL(ND,NP), the velocities.
%
% Input, real F(ND,NP), the forces.
%
% Input, real ACC(ND,NP), the accelerations.
%
% Input, real MASS, the mass of each particle.
%
% Input, real DT, the time step.
%
% Output, real POS(ND,NP), the updated positions.
%
% Output, real VEL(ND,NP), the updated velocities.
%
% Output, real ACC(ND,NP), the updated accelerations.
%
rmass = 1.0 / mass;

pos(1:nd,1:np) = pos(1:nd,1:np) + vel(1:nd,1:np) * dt ...
+ 0.5 * acc(1:nd,1:np) * dt * dt;

vel(1:nd,1:np) = vel(1:nd,1:np) ...
+ 0.5 * dt * ( f(1:nd,1:np) * rmass + acc(1:nd,1:np) );

acc(1:nd,1:np) = f(1:nd,1:np) * rmass;

return
end

Before you make changes to the source code as provid ed to you ensure th at t does what t is supposed to do. Once this is demonstrated, you need to do the following for a 3-D system with 250 particles 1. Run the simulation for 10000 time-steps. Pick an atom at random and record its displacement at every time step. Then carry out a Fourier transform of the atom's displacement in order to obtain the characteristic frequencies of vibration associated with the atom motion. (to be 2. submitted by (11/30) periodic boundary conditions for a simple square lattice of size 16x16 mplement a code with particles 3. Modify the interatomic potential to a harmonic potential, where each particle exclusively interacts with its first-nearest neighbors; use a spring constant of one. Provide 'small random displacements to each particle as an initial condition and repeat Step 1. (results to be presented in class on (12/5 and 12/7)

Explanation / Answer


np = 1000;
step_num = 500;
nd = 3;
dt = 0.0001;
mass = 1.0;

fprintf ( 1, ' ' );
fprintf ( 1, 'MD_PARALLEL ' );
fprintf ( 1, ' MATLAB version ' );
fprintf ( 1, ' ' );
fprintf ( 1, ' A molecular dynamics program. ' );
fprintf ( 1, ' ' );
fprintf ( 1, ' NP, the number of particles is %d. ', np );
fprintf ( 1, ' STEP_NUM, the number of time steps, is %d. ', step_num );
fprintf ( 1, ' DT, the time step size, is %f seconds. ', dt );
%
% Set the dimensions of the box.
%
box(1:nd) = 10.0;
%
% Set initial positions, velocities, and accelerations.
%
fprintf ( 1, ' ' );
fprintf ( 1, ' Initializing positions, velocities, and accelerations. ' );

[ pos, vel, acc ] = initialize ( np, nd, box );
%
% Compute the forces and energies.
%
fprintf ( 1, ' ' );
fprintf ( 1, ' Computing initial forces and energies. ' );

[ force, potential, kinetic ] = compute ( np, nd, pos, vel, mass );
%
% Save the initial total energy for use in the accuracy check.
%
e0 = potential + kinetic;
%
% This is the main time stepping loop:
% Compute forces and energies,
% Update positions, velocities, accelerations.
%
fprintf ( 1, ' ' );
fprintf ( 1, ' At each step, we report the potential and kinetic energies. ' );
fprintf ( 1, ' The sum of these energies should be a constant. ' );
fprintf ( 1, ' As an accuracy check, we also print the relative error ' );
fprintf ( 1, ' in the total energy. ' );
fprintf ( 1, ' ' );
fprintf ( 1, ' Step Potential Kinetic (P+K-E0)/E0 ' );
fprintf ( 1, ' Energy P Energy K Relative Energy Error ' );
fprintf ( 1, ' ' );
  
step_print_index = 0;
step_print_num = 10;
  
step = 0;
fprintf ( 1, ' %8d %14f %14f %14e ', ...
step, potential, kinetic, ( potential + kinetic - e0 ) / e0 );
step_print_index = step_print_index + 1;
step_print = floor ( ( step_print_index / step_print_num ) * step_num );
  
for step = 1 : step_num

[ force, potential, kinetic ] = compute ( np, nd, pos, vel, mass );

if ( step == step_print )
energy_lost = ( potential + kinetic - e0 ) / e0 ;
fprintf ( 1, ' %8d %14f %14f %14e ', ...
step, potential, kinetic, ( potential + kinetic - e0 ) / e0 );
step_print_index = step_print_index + 1;
step_print = floor ( ( step_print_index / step_print_num ) * step_num );
end
  
[ pos, vel, acc ] = update ( np, nd, pos, vel, force, acc, mass, dt );

end

fprintf ( 1, ' ' );
fprintf ( 1, 'MD_PARALLEL ' );
fprintf ( 1, ' Normal end of execution. ' );


return
end
function [ f, pot, kin ] = compute ( np, nd, pos, vel, mass )

=======================================================================================
%
% 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, PI2 ) ) )^2
%
% The derivative of the potential is:
%
% dv(x) = 2.0 * sin ( min ( x, PI2 ) ) * cos ( min ( x, PI2 ) )
% = sin ( 2.0 * min ( x, PI2 ) )
%

f = zeros ( nd, np );

pot = 0.0;

pi2 = pi / 2.0;

parfor i = 1 : np

Ri = pos - repmat ( pos( :, i ), 1, np ); % array of vectors to 'i'

D = sqrt ( sum ( Ri.^2 ) ); % array of distances

Ri = Ri( :, ( D > 0.0 ) );

D = D( D > 0.0 ); % save only pos values

D2 = D .* ( D <= pi2 ) + pi2 * ( D > pi2 ); % truncate the potential.

pot = pot + 0.5 * sum ( sin ( D2 ).^2 ); % accumulate pot. energy

f( :, i) = Ri * ( sin ( 2*D2 ) ./ D )'; % force on particle 'i'

end
%
% Compute kinetic energy.
%
kin = 0.5 * mass * sum ( diag ( vel' * vel ) );

return
end
function [ pos, vel, acc, seed ] = initialize ( np, nd, box, seed )

===================================================================================
for i = 1 : nd
pos(i,1:np) = box(i) * pos(i,1:np);
end
%
% Velocities and accelerations begin at 0.
%
vel(1:nd,1:np) = 0.0;
acc(1:nd,1:np) = 0.0;

return
end
function [ pos, vel, acc ] = update ( np, nd, pos, vel, f, acc, mass, dt )

========================================================================
% x(t+dt) = x(t) + v(t) * dt + 0.5 * a(t) * dt * dt
% v(t+dt) = v(t) + 0.5 * ( a(t) + a(t+dt) ) * dt
% a(t+dt) = f(t) / m
%


rmass = 1.0 / mass;

pos(1:nd,1:np) = pos(1:nd,1:np) + vel(1:nd,1:np) * dt ...
+ 0.5 * acc(1:nd,1:np) * dt * dt;

vel(1:nd,1:np) = vel(1:nd,1:np) ...
+ 0.5 * dt * ( f(1:nd,1:np) * rmass + acc(1:nd,1:np) );

acc(1:nd,1:np) = f(1:nd,1:np) * rmass;

return
end

Related Questions

full blank with correct answer US Airways Group 8688.0 International Paper 19500
Question #3258397
full credit far nos-mutiple chaice probienss.Follow the direction in the Witten
Question #3313491
full data please help this is my last lab ! (enthalpy change for the decompositi
Question #872578
full detail 17 The Bedouin will old Bedouin died and left 17 camels to his three
Question #3423565
full documentation here; http://cs.boisestate.edu/~cs121/projects/p4/ Here is my
Question #3731802
full employment level of GDP Which of the following is a mechanism Keynes advoca
Question #1215754
full points to anyone who answers in detail all of the questions An infinitely l
Question #2112970
full points to anyone who answers in detail all of the questions Two identically
Question #2112971

Navigate

Previous
function md ( nd, np, step_num, dt ) %******************************************
Next
function md ( nd, np, step_num, dt ) %******************************************

Integrity-first tutoring: explanations and feedback only — we do not complete graded work. Learn more.