Hi,
Just joined these forums as I've began to model materials engineering situations in C as part of my course.
Essentially I need to edit this program below so that it will produce a movie of the file which will then be played in JMol.
Here's the program as it currently stands:
What I need to do is modify the computer code so that it writes out the average position of the atoms in each block of size NBIN after each time interval of length TDUMP.Code:#include <stdio.h> #include <math.h> /* Simulation parameters */ #define N 5000 /* Total number of atoms: projectile + armour */ #define NP 500 /* Size of projectile */ #define V0 1.0e-2 /* Initial velocity of projectile in A/fs */ #define TSIM 2.0e+5 /* Length of simulation in fs */ #define DT 1.0 /* Time step in fs */ #define NBIN 100 /* Size of bin for output */ #define TDUMP 1000.0 /* Interval between dumps to file in fs */ /* Scales */ #define MASS 55.845 /* Mass in amu */ #define R0 1.4333 /* Bond length in A - fixed by (100) interplanar spacing */ #define E0 4.29 /* Bond strength in eV - from cohesive energy */ #define L0 1.12 /* Potential decay rate in A - from speed of sound */ #define VOLUME 5.8884 /* Volume per atom in A^3 */ int main () { int i, j, k, n, nsteps, ndump; double t, mass, aij, fac, fac2, PE, KE, U; double xp, vi, vbin, kebin, stressbin; double xm[N], x[N], a[N]; FILE *vcmfp, *tempfp, *stressfp; /* Convert mass into consistent units: [mass] = [E] * [t]^2 / [l]^2 = eV fs^2 / A^2 = 1.602176e-19 * 1.0e-15^2 / (1.0e-10^2 * 1.660538e-27) amu = 9.6485e-3 amu */ mass = MASS / 9.6485e-3; /* Initial conditions */ for (i = 0; i < N; i++) { x[i] = ((double) i) * R0; if (i < NP) xm[i] = x[i] - V0 * DT; else xm[i] = x[i]; a[i] = 0.0; } /* Get number of MD steps */ nsteps = (int)(TSIM/DT) + 1; ndump = (int)(TDUMP/DT); /* Open output files */ vcmfp = fopen ("velocity-cm.dat", "w"); tempfp = fopen ("temperature.dat", "w"); stressfp = fopen ("stress.dat", "w"); /* Main loop Use Morse potential: v(x) = E0 (exp(-2(x-R0)/L0) - 2 exp(-(x-R0)/L0)) F(x) = -(2 E0/L0) (exp(-2(x-R0)/L0) - exp(-(x-R0)/L0)) */ t = 0.0; for (n = 0; n < nsteps; n++) { /* Update output */ if (n%ndump == 0) { /* Trace */ PE = 0.0; vi = (x[N-1] - xm[N-1])/DT + 0.5*a[N-1]*DT; KE = 0.5 * mass * vi*vi; for (i = 0; i < N-1; i++) { fac = exp(-(x[i+1]-x[i]-R0)/L0); fac2 = fac * fac; PE += E0 * (fac2 - 2.0 * fac); vi = (x[i] - xm[i])/DT + 0.5*a[i]*DT; KE += 0.5 * mass * vi*vi; } U = PE + KE; printf ("%6.2f%% Energy = %10.2f\n", (t/TSIM)*100.0, U); /* Temperature and velocity of the centre of mass */ for (i = 0; i < N; i += NBIN) { vbin = 0.0; kebin = 0.0; for (j = 0; j < NBIN; j++) { k = i+j; vi = (x[k] - xm[k])/DT + 0.5*a[k]*DT; vbin += vi; kebin += 0.5 * mass * vi*vi; } fprintf (vcmfp, "%7.1f %6i %8.1f\n", t, i+1, vbin * 1.0e+5 / ((double) NBIN)); /* Velocity of CM in m/s */ kebin -= 0.5 * mass * vbin*vbin / ((double) NBIN); fprintf (tempfp, "%7.1f %6i %8.1f\n", t, i+1, kebin * 11604.0 / ((double) NBIN)); /* Temperature in K */ } fprintf (vcmfp, "\n"); fprintf (tempfp, "\n"); /* Stress */ for (i = 0; i < N; i += NBIN) { stressbin = 0.0; for (j = 0; j < NBIN; j++) { k = i+j; if (k < N-1) { fac = exp(-(x[k+1]-x[k]-R0)/L0); fac2 = fac * fac; stressbin += -(x[k+1]-x[k]) * (2.0*E0/L0) * (fac2 - fac); } } fprintf (stressfp, "%7.1f %6i %8.1f\n", t, i+1, stressbin * 160.218 / (VOLUME * (double) NBIN)); /* Stress in GPa */ } fprintf (stressfp, "\n"); } /* Get acceleration at time t */ a[0] = 0.0; for (i = 0; i < N-1; i++) { /* Morse factors */ fac = exp(-(x[i+1]-x[i]-R0)/L0); fac2 = fac * fac; /* Update force for bond */ aij = -(2.0 * E0 / (L0 * mass)) * (fac2 - fac); a[i] += aij; a[i+1] = -aij; } /* Advance positions and velocities to time t+dt using Verlet: x(t+dt) = 2 x(t) - x(t-dt) + a(t) dt^2 */ for (i = 0; i < N; i++) { xp = 2.0 * x[i] - xm[i] + a[i] * DT * DT; xm[i] = x[i]; x[i] = xp; } t += DT; } /* Close output files */ fclose (vcmfp); fclose (tempfp); fclose (stressfp); return 0; }
From what I'm assuming is that I'll need three loops, the first loop doing the timestep iterated by TDUMP, the next loop doing the set, having to add NBIN each time so that i knows where to start and the final loop doing the average. However I'm not really sure how to put this together.
I was looking at another program which created a video and have been trying to understand the steps in this one to apply it to the one above:
Any suggestions would be great,Code:#include <stdio.h> #include <stdlib.h> #include <math.h> #define PI 3.14159265358 int main (int *argc, char **argv) { int i, natom; double t, t_min, t_max, dt, period; double *x0, xi, amplitude, wavelength, bondlength; FILE *fp; /* Build chain of atoms */ natom = 10; bondlength = 2.0; x0 = (double *) malloc (natom * sizeof (double)); x0[0] = 0.0; for (i = 1; i < natom; i++) x0[i] = x0[i-1] + bondlength; /* Define oscillation parameters */ amplitude = 0.2; period = 100.0; wavelength = 5.0 * bondlength; /* Open output file */ fp = fopen ("movie.xyz", "w"); /* Write out atomic positions for each time */ t_min = 0.0; t_max = 5.0 * period; dt = (t_max - t_min) / 100.0; for (t = t_min; t <= t_max; t += dt) { fprintf (fp, "%i\n", natom); fprintf (fp, "t = %15.7e\n", t); for (i = 0; i < natom; i++) { xi = x0[i] + amplitude * sin(2.0*PI * (t / period - x0[i] / wavelength)); fprintf (fp, "Fe %15.7e %15.7e %15.7e\n", xi, 0.0, 0.0); } } /* Close output file and free memory */ fclose (fp); free (x0); return 0; }
thanks very much.
By the way I'm using Visual Studio to compile as that is on the computers in College.
