Having problems generating a video (.xyz)

This is a discussion on Having problems generating a video (.xyz) within the C Programming forums, part of the General Programming Boards category; Hi, Just joined these forums as I've began to model materials engineering situations in C as part of my course. ...

  1. #1
    Registered User
    Join Date
    Feb 2009
    Posts
    3

    Having problems generating a video (.xyz)

    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:

    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;
    }
    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.

    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:

    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;
    }
    Any suggestions would be great,

    thanks very much.

    By the way I'm using Visual Studio to compile as that is on the computers in College.

  2. #2
    Registered User
    Join Date
    Feb 2009
    Posts
    3
    Right, i've actually managed to do this myself in the end. However when building it the build succeeds but I get the old microsoft "it encountered a problem and needs to close".

    The new code is below and I was wondering if anyone could see if there is anything that could cause this, thanks.

    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 */
    #define NA	  50	   /* Number of atoms for movie file*/
    
    /* 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, xavp;
    	double xm[N], x[N], a[N];
    	FILE *vcmfp, *tempfp, *stressfp, *moviefp;
    
    	/* 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");
    	moviefp = fopen ("movie.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");
    
    		    /* moviefile */
    
    			fprintf (moviefp, "%i\n\n", NA);
    
    			for (i = 0; i < N; i += NBIN)
    			    {
    				xavp = 0;	
    				
    			/*Set ave. x position to zero*/
    				
    			for (j = 0; j < NBIN; j++)
    				{
    					k = i+k;
    					xavp += x[k];
    				}
    				xavp /= (double)(NBIN * NBIN);
    
    			/*Produce average x value then divide by NBIN to reduce the spacing*/
    
    				fprintf (moviefp, "Fe %f 0.00 0.00\n", xavp);
    			}
    		}
    
    		/* 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);
    	fclose (moviefp);
    	
    
    	return 0;
    }

  3. #3
    and the Hat of Guessing tabstop's Avatar
    Join Date
    Nov 2007
    Posts
    14,185
    I compiled your code and ran it in a debugger. Your program died on line 138, with a value of k equal to 10299 (which is certainly an out-of-bounds access). The line before that is k = i+k, which is different from all those other lines of k = i+j.

  4. #4
    Registered User
    Join Date
    Feb 2009
    Posts
    3
    Thanks very much!

Popular pages Recent additions subscribe to a feed

Similar Threads

  1. Doxygen failing
    By Elysia in forum A Brief History of Cprogramming.com
    Replies: 19
    Last Post: 04-16-2008, 01:24 PM
  2. Violent video games?
    By VirtualAce in forum A Brief History of Cprogramming.com
    Replies: 58
    Last Post: 04-26-2006, 01:43 PM
  3. Problem With My Box
    By HaVoX in forum Tech Board
    Replies: 9
    Last Post: 10-15-2005, 07:38 AM
  4. Quick Video Card Question
    By XSquared in forum Tech Board
    Replies: 3
    Last Post: 09-13-2005, 08:50 PM

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21