[Swift-commit] r8124 - SwiftTutorials/ATPESC_2014-08-14/swift-t/examples/05-md
wozniak at ci.uchicago.edu
wozniak at ci.uchicago.edu
Wed Aug 13 10:46:03 CDT 2014
Author: wozniak
Date: 2014-08-13 10:46:03 -0500 (Wed, 13 Aug 2014)
New Revision: 8124
Added:
SwiftTutorials/ATPESC_2014-08-14/swift-t/examples/05-md/md.c
Log:
Copied from:
URL: https://svn.ci.uchicago.edu/svn/vdl2/demo/xsede.2014.0425/src/md
Revision: 8122
Added: SwiftTutorials/ATPESC_2014-08-14/swift-t/examples/05-md/md.c
===================================================================
--- SwiftTutorials/ATPESC_2014-08-14/swift-t/examples/05-md/md.c (rev 0)
+++ SwiftTutorials/ATPESC_2014-08-14/swift-t/examples/05-md/md.c 2014-08-13 15:46:03 UTC (rev 8124)
@@ -0,0 +1,787 @@
+# include <stdlib.h>
+# include <stdio.h>
+# include <time.h>
+# include <math.h>
+
+int main ( int argc, char *argv[] );
+void compute ( int np, int nd, double pos[], double vel[],
+ double mass, double f[], double *pot, double *kin );
+double cpu_time ( void );
+double dist ( int nd, double r1[], double r2[], double dr[] );
+void initialize ( int np, int nd, double box[], int *seed, double pos[],
+ double vel[], double acc[] );
+double r8_uniform_01 ( int *seed );
+void timestamp ( void );
+void update ( int np, int nd, double pos[], double vel[], double f[],
+ double acc[], double mass, double dt );
+void snap ( int np, int nd, double pos[], double vel[], double f[],
+ double acc[], double mass, double dt );
+
+double scale_factor = 2.5, scale_offset = -2.0;
+char *printinfo = "0.05 1.0 0.2 0.05 50.0 0.1";
+char *outfile = "md.dat";
+char *trjfile = "md.trj.tgz";
+
+
+/******************************************************************************/
+
+int main ( int argc, char *argv[] )
+
+/******************************************************************************/
+/*
+ Purpose:
+
+ 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.
+
+ Usage:
+
+ md nd np step_num print_step_num dt mass printinfo scale_factor scale_offset seed outFile trajectoryFile
+ 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);
+ * print_step_num is the number of snapshot prints (10 for instance);
+ * dt is size of timestep;
+ * mass is particle mass;
+ * printinfo is a string to append to each particle coord
+ * scale_offset and scale_factor are used to scale particle positions for logging/rendering (FIXME)
+ * seed sets the initial configuration
+
+
+ Licensing:
+
+ This code is distributed under the GNU LGPL license.
+
+ Modified:
+
+ 05 November 2010
+
+ Author:
+
+ Original FORTRAN90 version by Bill Magro.
+ C version by John Burkardt.
+
+ Parameters:
+
+ None
+*/
+{
+ double *acc;
+ double *box;
+ double ctime;
+ double ctime1;
+ double ctime2;
+ double dt = 0.0001;
+ double e0;
+ double *force;
+ int i;
+ int id;
+ double kinetic;
+ double mass = 1.0 * .0001;
+ int nd;
+ int np;
+ double *pos;
+ double potential;
+ int seed = 123456789;
+ int step;
+ int step_num;
+ int step_print;
+ int step_print_index;
+ int step_print_num=10;
+ double *vel;
+
+ timestamp ( );
+ printf ( "\n" );
+ printf ( "MD\n" );
+ printf ( " C version\n" );
+ printf ( " A molecular dynamics program.\n" );
+/*
+ Get the spatial dimension.
+*/
+ if ( 1 < argc )
+ {
+ nd = atoi ( argv[1] );
+ }
+ else
+ {
+ printf ( "\n" );
+ printf ( " Enter ND, the spatial dimension (2 or 3).\n" );
+ scanf ( "%d", &nd );
+ }
+//
+// Get the number of points.
+//
+ if ( 2 < argc )
+ {
+ np = atoi ( argv[2] );
+ }
+ else
+ {
+ printf ( "\n" );
+ printf ( " Enter NP, the number of points (500, for instance).\n" );
+ scanf ( "%d", &np );
+ }
+//
+// Get the number of time steps.
+//
+ if ( 3 < argc )
+ {
+ step_num = atoi ( argv[3] );
+ }
+ else
+ {
+ printf ( "\n" );
+ printf ( " Enter ND, the number of time steps (500 or 1000, for instance).\n" );
+ scanf ( "%d", &step_num );
+ }
+ /*
+ Get any additional args (command-line only)
+ md nd np step_num [ step__print_num dt mass printinfo scale_factor scale_offset randomseed outfile trjfile ]
+ */
+ if ( 4 < argc )
+ {
+ step_print_num = atoi ( argv[4] );
+ }
+ if ( 5 < argc )
+ {
+ dt = atof ( argv[5] );
+ }
+ if ( 6 < argc )
+ {
+ mass = atof ( argv[6] );
+ }
+ if ( 7 < argc )
+ {
+ printinfo = ( argv[7] );
+ }
+ if ( 8 < argc )
+ {
+ scale_factor = atof ( argv[8] );
+ }
+ if ( 9 < argc )
+ {
+ scale_offset = atof ( argv[9] );
+ }
+ if ( 10 < argc )
+ {
+ seed = atof ( argv[10] );
+ }
+ if ( 11 < argc )
+ {
+ outfile = argv[11];
+ }
+ if ( 12 < argc )
+ {
+ trjfile = argv[12];
+ }
+
+/*
+ Report.
+*/
+ printf ( "\n" );
+ printf ( " MD: Argument count: %d\n", argc );
+ printf ( " ND, the spatial dimension, is %d\n", nd );
+ printf ( " NP, the number of particles in the simulation, is %d\n", np );
+ printf ( " STEP_NUM, the number of time steps, is %d\n", step_num );
+ printf ( " STEP_PRINT_NUM, the number of snapshots to print, is %d\n", step_print_num );
+ printf ( " DT, the size of each time step, is %f\n", dt );
+ printf ( " MASS, the particle mass, is %f\n", mass );
+ printf ( " PRINTINFO, the pass-through info to c-ray, is %s\n", printinfo );
+ printf ( " SCALE_FACTOR, the particle position scaling factor, is %f\n", scale_factor );
+ printf ( " SCALE_OFFSET, the particle position scaling offset, is %f\n", scale_offset );
+ printf ( " SEED, the simulation randomization seed, is %d\n", seed );
+/*
+ Allocate memory.
+*/
+ acc = ( double * ) malloc ( nd * np * sizeof ( double ) );
+ box = ( double * ) malloc ( nd * sizeof ( double ) );
+ force = ( double * ) malloc ( nd * np * sizeof ( double ) );
+ pos = ( double * ) malloc ( nd * np * sizeof ( double ) );
+ vel = ( double * ) malloc ( nd * np * sizeof ( double ) );
+/*
+ Set the dimensions of the box.
+*/
+ for ( i = 0; i < nd; i++ )
+ {
+ box[i] = 10.0;
+ }
+
+ printf ( "\n" );
+ printf ( " Initializing positions, velocities, and accelerations.\n" );
+/*
+ Set initial positions, velocities, and accelerations.
+*/
+ initialize ( np, nd, box, &seed, pos, vel, acc );
+/*
+ Compute the forces and energies.
+*/
+ printf ( "\n" );
+ printf ( " Computing initial forces and energies.\n" );
+
+ compute ( np, nd, pos, vel, mass, force, &potential, &kinetic );
+
+ e0 = potential + kinetic;
+/*
+ This is the main time stepping loop:
+ Compute forces and energies,
+ Update positions, velocities, accelerations.
+*/
+ printf ( "\n" );
+ printf ( " At each step, we report the potential and kinetic energies.\n" );
+ printf ( " The sum of these energies should be a constant.\n" );
+ printf ( " As an accuracy check, we also print the relative error\n" );
+ printf ( " in the total energy.\n" );
+ printf ( "\n" );
+ printf ( " Step Potential Kinetic (P+K-E0)/E0\n" );
+ printf ( " Energy P Energy K Relative Energy Error\n" );
+ printf ( "\n" );
+
+ FILE *ofile = fopen(outfile,"w");
+ fprintf (ofile, " Step Potential Kinetic RelativeErr\n" );
+
+ step_print = 0;
+ step_print_index = 0;
+
+ step = 0;
+ printf ( " %8d %14f %14f %14e\n",
+ step, potential, kinetic, ( potential + kinetic - e0 ) / e0 );
+ fprintf ( ofile, " %8d %14f %14f %14e\n",
+ step, potential, kinetic, ( potential + kinetic - e0 ) / e0 );
+ step_print_index = step_print_index + 1;
+ step_print = ( step_print_index * step_num ) / step_print_num;
+
+ ctime1 = cpu_time ( );
+
+ for ( step = 1; step <= step_num; step++ )
+ {
+ compute ( np, nd, pos, vel, mass, force, &potential, &kinetic );
+
+ if ( step == step_print )
+ {
+ printf ( " %8d %14f %14f %14e\n", step, potential, kinetic,
+ ( potential + kinetic - e0 ) / e0 );
+ fprintf ( ofile, " %8d %14f %14f %14e\n", step, potential, kinetic,
+ ( potential + kinetic - e0 ) / e0 );
+ step_print_index = step_print_index + 1;
+ step_print = ( step_print_index * step_num ) / step_print_num;
+ snap ( np, nd, pos, vel, force, acc, mass, dt );
+ }
+ update ( np, nd, pos, vel, force, acc, mass, dt );
+ }
+ ctime2 = cpu_time ( );
+ ctime = ctime2 - ctime1;
+
+ printf ( "\n" );
+ printf ( " Elapsed cpu time for main computation:\n" );
+ printf ( " %f seconds.\n", ctime );
+
+ free ( acc );
+ free ( box );
+ free ( force );
+ free ( pos );
+ free ( vel );
+ char tarcmd[2000];
+ sprintf(tarcmd,"tar zcf %s md??.trj",trjfile);
+ system(tarcmd);
+/*
+ Terminate.
+*/
+ printf ( "\n" );
+ printf ( "MD\n" );
+ printf ( " Normal end of execution.\n" );
+
+ printf ( "\n" );
+ timestamp ( );
+
+ fclose(ofile);
+ return 0;
+}
+/******************************************************************************/
+
+void compute ( int np, int nd, double pos[], double vel[],
+ double mass, double f[], double *pot, double *kin )
+
+/******************************************************************************/
+/*
+ Purpose:
+
+ 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, 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 ) )
+
+ Licensing:
+
+ This code is distributed under the GNU LGPL license.
+
+ Modified:
+
+ 21 November 2007
+
+ Author:
+
+ Original FORTRAN90 version by Bill Magro.
+ C version by John Burkardt.
+
+ Parameters:
+
+ Input, int NP, the number of particles.
+
+ Input, int ND, the number of spatial dimensions.
+
+ Input, double POS[ND*NP], the position of each particle.
+
+ Input, double VEL[ND*NP], the velocity of each particle.
+
+ Input, double MASS, the mass of each particle.
+
+ Output, double F[ND*NP], the forces.
+
+ Output, double *POT, the total potential energy.
+
+ Output, double *KIN, the total kinetic energy.
+*/
+{
+ double d;
+ double d2;
+ int i;
+ int j;
+ int k;
+ double ke;
+ double pe;
+ double PI2 = 3.141592653589793 / 2.0;
+ double rij[3];
+
+ pe = 0.0;
+ ke = 0.0;
+
+ for ( k = 0; k < np; k++ )
+ {
+/*
+ Compute the potential energy and forces.
+*/
+ for ( i = 0; i < nd; i++ )
+ {
+ f[i+k*nd] = 0.0;
+ }
+
+ for ( j = 0; j < np; j++ )
+ {
+ if ( k != j )
+ {
+ d = dist ( nd, pos+k*nd, pos+j*nd, rij );
+/*
+ Attribute half of the potential energy to particle J.
+*/
+ if ( d < PI2 )
+ {
+ d2 = d;
+ }
+ else
+ {
+ d2 = PI2;
+ }
+
+ pe = pe + 0.5 * pow ( sin ( d2 ), 2 );
+
+ for ( i = 0; i < nd; i++ )
+ {
+ f[i+k*nd] = f[i+k*nd] - rij[i] * sin ( 2.0 * d2 ) / d;
+ }
+ }
+ }
+/*
+ Compute the kinetic energy.
+*/
+ for ( i = 0; i < nd; i++ )
+ {
+ ke = ke + vel[i+k*nd] * vel[i+k*nd];
+ }
+ }
+
+ ke = ke * 0.5 * mass;
+
+ *pot = pe;
+ *kin = ke;
+
+ return;
+}
+/*******************************************************************************/
+
+double cpu_time ( void )
+
+/*******************************************************************************/
+/*
+ Purpose:
+
+ CPU_TIME reports the total CPU time for a program.
+
+ Licensing:
+
+ This code is distributed under the GNU LGPL license.
+
+ Modified:
+
+ 27 September 2005
+
+ Author:
+
+ John Burkardt
+
+ Parameters:
+
+ Output, double CPU_TIME, the current total elapsed CPU time in second.
+*/
+{
+ double value;
+
+ value = ( double ) clock ( ) / ( double ) CLOCKS_PER_SEC;
+
+ return value;
+}
+/******************************************************************************/
+
+double dist ( int nd, double r1[], double r2[], double dr[] )
+
+/******************************************************************************/
+/*
+ Purpose:
+
+ DIST computes the displacement (and its norm) between two particles.
+
+ Licensing:
+
+ This code is distributed under the GNU LGPL license.
+
+ Modified:
+
+ 21 November 2007
+
+ Author:
+
+ Original FORTRAN90 version by Bill Magro.
+ C version by John Burkardt.
+
+ Parameters:
+
+ Input, int ND, the number of spatial dimensions.
+
+ Input, double R1[ND], R2[ND], the positions of the particles.
+
+ Output, double DR[ND], the displacement vector.
+
+ Output, double D, the Euclidean norm of the displacement.
+*/
+{
+ double d;
+ int i;
+
+ d = 0.0;
+ for ( i = 0; i < nd; i++ )
+ {
+ dr[i] = r1[i] - r2[i];
+ d = d + dr[i] * dr[i];
+ }
+ d = sqrt ( d );
+
+ return d;
+}
+/******************************************************************************/
+
+void initialize ( int np, int nd, double box[], int *seed, double pos[],
+ double vel[], double acc[] )
+
+/******************************************************************************/
+/*
+ Purpose:
+
+ INITIALIZE initializes the positions, velocities, and accelerations.
+
+ Licensing:
+
+ This code is distributed under the GNU LGPL license.
+
+ Modified:
+
+ 20 July 2008
+
+ Author:
+
+ Original FORTRAN90 version by Bill Magro.
+ C version by John Burkardt.
+
+ Parameters:
+
+ Input, int NP, the number of particles.
+
+ Input, int ND, the number of spatial dimensions.
+
+ Input, double BOX[ND], specifies the maximum position
+ of particles in each dimension.
+
+ Input, int *SEED, a seed for the random number generator.
+
+ Output, double POS[ND*NP], the position of each particle.
+
+ Output, double VEL[ND*NP], the velocity of each particle.
+
+ Output, double ACC[ND*NP], the acceleration of each particle.
+*/
+{
+ int i;
+ int j;
+/*
+ Give the particles random positions within the box.
+*/
+ for ( j = 0; j < np; j++ )
+ {
+ for ( i = 0; i < nd; i++ )
+ {
+ pos[i+j*nd] = box[i] * r8_uniform_01 ( seed );
+ vel[i+j*nd] = 0.0;
+ acc[i+j*nd] = 0.0;
+ }
+ }
+ return;
+}
+/******************************************************************************/
+
+double r8_uniform_01 ( int *seed )
+
+/******************************************************************************/
+/*
+ Purpose:
+
+ R8_UNIFORM_01 is a unit pseudorandom R8.
+
+ Discussion:
+
+ This routine implements the recursion
+
+ seed = 16807 * seed mod ( 2^31 - 1 )
+ unif = seed / ( 2^31 - 1 )
+
+ The integer arithmetic never requires more than 32 bits,
+ including a sign bit.
+
+ Licensing:
+
+ This code is distributed under the GNU LGPL license.
+
+ Modified:
+
+ 11 August 2004
+
+ Author:
+
+ John Burkardt
+
+ Reference:
+
+ Paul Bratley, Bennett Fox, Linus Schrage,
+ A Guide to Simulation,
+ Springer Verlag, pages 201-202, 1983.
+
+ Bennett Fox,
+ Algorithm 647:
+ Implementation and Relative Efficiency of Quasirandom
+ Sequence Generators,
+ ACM Transactions on Mathematical Software,
+ Volume 12, Number 4, pages 362-376, 1986.
+
+ Parameters:
+
+ Input/output, int *SEED, a seed for the random number generator.
+
+ Output, double R8_UNIFORM_01, a new pseudorandom variate, strictly between
+ 0 and 1.
+*/
+{
+ int k;
+ double r;
+
+ k = *seed / 127773;
+
+ *seed = 16807 * ( *seed - k * 127773 ) - k * 2836;
+
+ if ( *seed < 0 )
+ {
+ *seed = *seed + 2147483647;
+ }
+
+ r = ( double ) ( *seed ) * 4.656612875E-10;
+
+ return r;
+}
+/******************************************************************************/
+
+void timestamp ( void )
+
+/******************************************************************************/
+/*
+ Purpose:
+
+ TIMESTAMP prints the current YMDHMS date as a time stamp.
+
+ Example:
+
+ 31 May 2001 09:45:54 AM
+
+ Licensing:
+
+ This code is distributed under the GNU LGPL license.
+
+ Modified:
+
+ 24 September 2003
+
+ Author:
+
+ John Burkardt
+
+ Parameters:
+
+ None
+*/
+{
+# define TIME_SIZE 40
+
+ static char time_buffer[TIME_SIZE];
+ const struct tm *tm;
+ size_t len;
+ time_t now;
+
+ now = time ( NULL );
+ tm = localtime ( &now );
+
+ len = strftime ( time_buffer, TIME_SIZE, "%d %B %Y %I:%M:%S %p", tm );
+
+ printf ( "%s\n", time_buffer );
+
+ return;
+# undef TIME_SIZE
+}
+/******************************************************************************/
+
+void update ( int np, int nd, double pos[], double vel[], double f[],
+ double acc[], double mass, double dt )
+
+/******************************************************************************/
+/*
+ Purpose:
+
+ 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:
+
+ 21 November 2007
+
+ Author:
+
+ Original FORTRAN90 version by Bill Magro.
+ C version by John Burkardt.
+
+ Parameters:
+
+ Input, int NP, the number of particles.
+
+ Input, int ND, the number of spatial dimensions.
+
+ Input/output, double POS[ND*NP], the position of each particle.
+
+ Input/output, double VEL[ND*NP], the velocity of each particle.
+
+ Input, double F[ND*NP], the force on each particle.
+
+ Input/output, double ACC[ND*NP], the acceleration of each particle.
+
+ Input, double MASS, the mass of each particle.
+
+ Input, double DT, the time step.
+*/
+{
+ int i;
+ int j;
+ double rmass;
+
+ rmass = 1.0 / mass;
+
+ for ( j = 0; j < np; j++ )
+ {
+ for ( i = 0; i < nd; i++ )
+ {
+ pos[i+j*nd] = pos[i+j*nd] + vel[i+j*nd] * dt + 0.5 * acc[i+j*nd] * dt * dt;
+ vel[i+j*nd] = vel[i+j*nd] + 0.5 * dt * ( f[i+j*nd] * rmass + acc[i+j*nd] );
+ acc[i+j*nd] = f[i+j*nd] * rmass;
+ }
+ /*
+s -0.5 -0.3 -1 0.05 1.0 0.2 0.05 50.0 0.1
+s 0.5 -0.4 0 0.05 0.1 0.85 1.0 50.0 0.1
+ */
+ char *s = "0.05 1.0 0.2 0.05 50.0 0.1";
+ /* printf("s %f %f %f %s\n", pos[0+j*nd]/20.0, pos[1+j*nd]/20.0, pos[2+j*nd]/20.0, s); */
+ }
+
+ return;
+}
+
+
+double scale (double x)
+{
+ return( (x / scale_factor) - scale_offset);
+ /* return( (x / 2.5) - 2.0); */
+}
+
+static int snapid = 0;
+
+void snap ( int np, int nd, double pos[], double vel[], double f[],
+ double acc[], double mass, double dt )
+{
+ int j;
+
+ char snapfile[100];
+ sprintf(snapfile, "md%02d.trj", snapid);
+
+ FILE *sf = fopen(snapfile, "w");
+
+ for ( j = 0; j < np; j++ )
+ {
+ fprintf(sf, "s %f %f %f %s\n", scale(pos[0+j*nd]), scale(pos[1+j*nd]), scale(pos[2+j*nd]), printinfo);
+ }
+
+ fclose(sf);
+ snapid++;
+}
More information about the Swift-commit
mailing list