/*==========================================================================*/ /* */ /* Molecular Dynamics */ /* */ /* Constant Particle, Pressure, Energy */ /* */ /* (c) Dieter W. Heermann */ /* Heidelberg 2000 */ /* Version 31.05.2000 */ /* */ /*==========================================================================*/ /* ==== Table of constant values ==== */ # define N 256 /* Number of particles */ # define flag_velocity_scaling 1 /* Perform a scaling of */ /* the velocities */ /* ==== Include files for this program ==== */ # include # include # include # define RAND_MAX 2147483647 /*--------------------------------------------------------------------------*/ /* Main program */ /*--------------------------------------------------------------------------*/ int main(argc,argv) int argc; char *argv[]; { /*---------------------------------------------------------- */ /* Declarations */ /*---------------------------------------------------------- */ /* ==== First the arrays which we are using ==== */ float x[N], y[N], z[N]; float xdv[N], ydv[N], zdv[N]; float vhx[N], vhy[N], vhz[N]; float fx[N], fy[N], fz[N]; /* ==== Now the rest of the variables ==== */ int i,j,k,n,step, lg; /* loop indices */ int seed; int n_md; int f_return_code; int index; int stopScaling, analysis; int restart; float density; float side ; float sideh, rcoff, rcoffs; float h, hsq, hsq2; float px, py, pz; float tscale; float kinetic, potential, energy; float partialKinetic; float virial, pressure; float rp; float T1,tref; float temperature; float rd; float rd2,rd3,rd4,rd6,rd7; float xi, yi, zi; float kx, ky, kz; float r148; float tmp0,tmp1,tmp2; float tmp3,tmp4,tmp5; float lattice_const; float V, sideScale, pExt, massV; float vp,pd, vol; /* ==== Dcl's for the file system ==== */ FILE *config_new_cid; FILE *config_old_cid; /*---------------------------------------------------------- */ /* End of the Declarations */ /*---------------------------------------------------------- */ /*---------------------------------------------------------- */ /* Start with setting up procedures and perhaps */ /* read in the data from a previous run */ /*---------------------------------------------------------- */ /* ==== Print the header for this simulation ==== */ printf("\n\n"); printf("(c) Dieter W. Heermann \n\n"); printf("Version 0.8 (2000) of the program: md_npe \n\n"); printf("\n\n"); printf("------------------------------------------------------------ \n"); printf("Initial set-up procedure ... \n\n"); /* ==== Simulation parameters ==== */ tref = 2.53; pExt = 3.18; density = 0.636; massV = 0.01; rcoff = 2.5; h = 0.064; n_md = 200; seed = 1234; restart = 0; analysis = 1; stopScaling = 1000; /* ==== Prepare by initializing all necessary quantities ==== */ side = pow((double)N / density,(double)0.3333333); sideh = side * (float).5; sideScale = 1.0 / side; lattice_const = side / 4.0; V = side * side * side; hsq = h * h; hsq2 = hsq * (float).5; rcoffs = rcoff * rcoff; tscale = 16.0 / (N - 1.0); T1 = tref; printf("Side %f",side); printf("Opening file: ... "); config_new_cid = fopen("config_new.cfg","w"); if ( config_new_cid == NULL ) {printf("failed to open file \n"); exit(20); } else { printf("open \n"); } if ( restart == 0 ) { /* ==== Set up the initial configuration for the ==== */ /* ==== particles inside a cubical box ==== */ n = 0; for (lg = 0; lg < 2; lg++) { for (i=0; i < 4; i++) { for (j=0; j < 4; j++) { for (k=0; k < 4; k++) { x[n] = i * lattice_const + lg * lattice_const * 0.5; y[n] = j * lattice_const + lg * lattice_const * 0.5; z[n] = k * lattice_const; n++; } } } } for (lg = 1; lg < 3; lg++) { for (i=0; i < 4; i++) { for (j=0; j < 4; j++) { for (k=0; k < 4; k++) { x[n] = i * lattice_const + (2-lg) * lattice_const * 0.5; y[n] = j * lattice_const + (lg-1) * lattice_const * 0.5; z[n] = k * lattice_const + lattice_const * 0.5; n++; } } } } /* ==== Start with fresh velocities ==== */ for ( n=0; n < N; n++) { vhx[n] = ((float) rand() / RAND_MAX ) * ( 0.5 -((float)rand() / RAND_MAX )); vhy[n] = ((float) rand() / RAND_MAX ) * ( 0.5 -((float)rand() / RAND_MAX )); vhz[n] = ((float) rand() / RAND_MAX ) * ( 0.5 -((float)rand() / RAND_MAX )); } px = py = pz = 0.0; for (n=0; n < N; n++) { px += vhx[n]; py += vhy[n]; pz += vhz[n]; } px /= N; py /= N; pz /= N; for (n=0; n < N; n++) { vhx[n] -= px; vhy[n] -= py; vhz[n] -= pz; } kinetic = (float)0.; for (n = 0; n < N; n++) { kinetic += vhx[n] * vhx[n] + vhy[n] * vhy[n] + vhz[n] * vhz[n]; } kinetic /= hsq; for (n = 0; n < N; n++) { vhx[n] = (vhx[n] * sqrt(tref / (tscale * kinetic))) * sideScale; vhy[n] = (vhy[n] * sqrt(tref / (tscale * kinetic))) * sideScale; vhz[n] = (vhz[n] * sqrt(tref / (tscale * kinetic))) * sideScale; } } else { /* ==== Read in the previous configuration ==== */ printf("Opening file: ... "); config_old_cid = fopen("config_old.cfg","r"); if ( config_old_cid == NULL ) { printf("failed to open file \n"); exit(11); } else { printf("open \n"); } for(n=0; n < N; n++) { fscanf(config_old_cid,"%d %f %f %f %f %f %f\n", &index, &tmp0,&tmp1,&tmp2, &tmp3,&tmp4,&tmp5); x[n] = tmp0; y[n] = tmp1; z[n] = tmp2; vhx[n] = tmp3 * sideScale; vhy[n] = tmp4 * sideScale; vhz[n] = tmp5 * sideScale; } printf(" done\n"); f_return_code = fclose(config_old_cid); printf("have closed the file: return code %5d \n", f_return_code); } /* --- Config generation --- */ printf("done\n"); /* ==== Scale the absolute coordinates ==== */ for (n=0; n < N; n++) { xdv[n] = x[n] * sideScale; ydv[n] = y[n] * sideScale; zdv[n] = z[n] * sideScale; } /* ==== Compute the initial forces on the particles ==== */ printf("initialize forces ..."); potential = 0.0; virial = 0.0; virial = 0.0; for (n = 0; n < N; n++) { fx[n] = (float)0.; fy[n] = (float)0.; fz[n] = (float)0.; } /* ---- Compute the Lennard-Jones forces ---- */ for (i = 0; i < N; i++) { for (j = i+1; j < N; j++) { xi = x[i] - x[j]; yi = y[i] - y[j]; zi = z[i] - z[j]; /* ---- Minimum Image Convention ---- */ if (xi < -sideh) { xi += side;} if (xi > sideh) { xi -= side;} if (yi < -sideh) { yi += side;} if (yi > sideh) { yi -= side;} if (zi < -sideh) { zi += side;} if (zi > sideh) { zi -= side;} rd = xi*xi + yi*yi + zi*zi; if (rd < rcoffs) { rd2 = rd * rd; rd3 = rd * rd2; rd4 = rd2 * rd2; rd6 = rd2 * rd4; rd7 = rd3 * rd4; potential += ((1.0 /rd6) - (1.0 / rd3)); r148 = (1.0 / rd7) - (1.0 / rd4) * (float).5; virial -= rd * r148; kx = xi * r148; fx[i] += kx; fx[j] -= kx; ky = yi * r148; fy[i] += ky; fy[j] -= ky; kz = zi * r148; fz[i] += kz; fz[j] -= kz; } } } /* ==== Scale the forces properly ==== */ for (n = 0; n < N; n++) { fx[n] *= hsq2; fy[n] *= hsq2; fz[n] *= hsq2; } vp = 0.0; pd = 0.0; vol = 0.0; printf(" done \n\n"); printf("potential energy after set up: %f\n",potential); printf("done\n"); printf("finished with the preparations\n\n"); printf("----------------------------------------------------------\n\n"); /*---------------------------------------------------------- */ /* S i m u l a t i o n P a r t */ /*---------------------------------------------------------- */ for (step = 0; step < n_md; step++) { /* ==== Advance positions ==== */ for (n = 0; n < N; n++) { xdv[n] = xdv[n] + vhx[n] + fx[n] - vol * vhx[n]; ydv[n] = ydv[n] + vhy[n] + fy[n] - vol * vhy[n]; zdv[n] = zdv[n] + vhz[n] + fz[n] - vol * vhz[n]; } V += vp + pd; /* ==== Apply periodic boundary conditions ==== */ for (n = 0; n < N; n++) { if ( xdv[n] < 0 ) xdv[n] += 1.0; if ( xdv[n] > 1.0 ) xdv[n] -= 1.0; if ( ydv[n] < 0 ) ydv[n] += 1.0; if ( ydv[n] > 1.0 ) ydv[n] -= 1.0; if ( zdv[n] < 0 ) zdv[n] += 1.0; if ( zdv[n] > 1.0 ) zdv[n] -= 1.0; } /* ==== Compute the partial velocities ==== */ for (n = 0; n < N; n++) { vhx[n] += fx[n] -vol * vhx[n]; vhy[n] += fy[n] -vol * vhy[n]; vhz[n] += fz[n] -vol * vhz[n]; } /* ==== Compute the partial volume velocities ==== */ vp += pd; /* ==== Unscale the coordinates ==== */ side = pow(V, (1.0 / 3.0)); sideh = 0.5 * side; sideScale = 1.0 / side; for (n = 0; n < N; n++) { x[n] = xdv[n] * side; y[n] = ydv[n] * side; z[n] = zdv[n] * side; } /* ==== Compute the forces on the particles ==== */ virial = (float)0.; potential = (float)0.; for (n = 0; n < N; n++) { fx[n] = (float)0.; fy[n] = (float)0.; fz[n] = (float)0.; } /* ==== Compute the Lennard-Jones forces ==== */ for (i = 0; i < N; i++) { for (j = i+1; j < N; j++) { xi = x[i] - x[j]; yi = y[i] - y[j]; zi = z[i] - z[j]; /* ---- Minimum Image Convention ---- */ if (xi < -sideh) { xi += side;} if (xi > sideh) { xi -= side;} if (yi < -sideh) { yi += side;} if (yi > sideh) { yi -= side;} if (zi < -sideh) { zi += side;} if (zi > sideh) { zi -= side;} rd = xi*xi + yi*yi + zi*zi; if (rd < rcoffs) { rd2 = rd * rd; rd3 = rd * rd2; rd4 = rd2 * rd2; rd6 = rd2 * rd4; rd7 = rd3 * rd4; potential += ((1.0 /rd6) - (1.0 / rd3)); r148 = (1.0 / rd7) - (1.0 / rd4) * (float).5; virial -= rd * r148; kx = xi * r148; fx[i] += kx; fx[j] -= kx; ky = yi * r148; fy[i] += ky; fy[j] -= ky; kz = zi * r148; fz[i] += kz; fz[j] -= kz; } } } /* ==== scale the forces properly ==== */ for (n = 0; n < N; n++) { fx[n] *= hsq2 * sideScale; fy[n] *= hsq2 * sideScale; fz[n] *= hsq2 * sideScale; } partialKinetic = 0.0; for (n = 0; n < N; n++) { partialKinetic += vhx[n] * vhx[n] + vhy[n] * vhy[n] + vhz[n] * vhz[n]; } partialKinetic /= hsq; partialKinetic = partialKinetic * side * side; rp = 16.0 * (partialKinetic - virial) / V; pd = 0.5 * hsq * (rp -pExt) / massV; vp += pd; vol = 0.5 * vp / V; /* ==== Compute the velocities ==== */ for (n = 0; n < N; n++) { vhx[n] += fx[n] - vol * vhx[n]; vhy[n] += fy[n] - vol * vhy[n]; vhz[n] += fz[n] - vol * vhz[n]; } /* ==== Compute the kinetic energy ==== */ kinetic = (float)0.; for (n = 0; n < N; n++) { kinetic += vhx[n] * vhx[n] + vhy[n] * vhy[n] + vhz[n] * vhz[n]; } kinetic /= hsq; kinetic *= side * side; /* ==== This part will only be used for the set-up ==== */ /* ==== and initial thermalization process ==== */ if ( flag_velocity_scaling == 1) { if ( step < stopScaling ) { for (n = 0; n < N; n++) { vhx[n] = (vhx[n] * sqrt(tref / (tscale * kinetic))); vhy[n] = (vhy[n] * sqrt(tref / (tscale * kinetic))); vhz[n] = (vhz[n] * sqrt(tref / (tscale * kinetic))); } } } /* end flag for velcocity scaling */ /* ------------------------------------------------------------- */ /* ---- analysis part ---- */ /* ------------------------------------------------------------- */ temperature = tscale * kinetic; pressure = density* 16.0 * (kinetic - virial) / N; potential = 4.0 * potential; kinetic = 24.0 * kinetic; energy = kinetic + potential; if ( ( step /analysis ) * analysis == step ) { printf("%d %f %f %f %f %f %f\n", step, kinetic, potential, energy, temperature, pressure, V); } } /* ---- sweep loop ---- */ /*-------------------------------------------------------------------*/ /* Finish up everything */ /*-------------------------------------------------------------------*/ printf("writing file: ... "); for(n=0; n < N; n++) { tmp0 = x[n]; tmp1 = y[n]; tmp2 = z[n]; tmp3 = vhx[n]; tmp4 = vhy[n]; tmp5 = vhz[n]; fprintf(config_new_cid," %d %f %f %f %f %f %f\n", n,tmp0,tmp1,tmp2,tmp3,tmp4,tmp5); } printf("done\n"); return 0; } /* ----------END OF MAIN PROGRAM ------------*/