strnsub.c

一个很好的分子动力学程序

#include <matrix.h>	/* ludcmp  lubksp  */
#include <stdio.h>	/* printf			 */

/*  LIBRARY FUNCTION PROTOTYPES  */
void exit (int);
double fabs (double);
double sqrt (double);

void inverse	  (float *qxx, float *qxy, float *qxz,
						float *qyy, float *qyz, float *qzz) {
	float a[3][3], y[3][3], d;
	int	ind[3], i,j;
	a[0][0] = *qxx;
	a[0][1] = *qxy;
	a[0][2] = *qxz;
	a[1][0] = *qxy;
	a[1][1] = *qyy;
	a[1][2] = *qyz;
	a[2][0] = *qxz;
	a[2][1] = *qyz;
	a[2][2] = *qzz;
	for (i=0;i<3;i++)
	{
		for (j=0;j<3;j++)
			y[i][j] = 0;
		y[i][i] = 1;
	}
	ludcmp ( (float *) a, 3, 3, ind, &d);
	for (j=0;j<3;j++)
		lubksb ( (float *) a, 3, 3, ind, y[j]);
	*qxx = y[0][0];
	*qxy = y[0][1];
	*qxz = y[0][2];
	*qyy = y[1][1];
	*qyz = y[1][2];
	*qzz = y[2][2];
}


float sign (float a, float b) {
	if (b<0)
		return (-fabs(a));
	else
		return(fabs(a));
}


void tred2 (float a[3][3], int n, float *d, float *e) {
	int i,j,k,l;
	float scale,hh,h,g,f;

	if (n>1) {
		for (i=n-1;i>0;i--) {
			l = i-1;
			h = 0.0;
			if (l>0) {
				scale = 0.0;
				for (k=0;k<=l;k++)
					scale += fabs(a[i][k]);
				if (scale==0.0)
					e[i] = a[i][l];
				else {
					for (k=0;k<=l;k++)  {
						a[i][k] /= scale;
						h += a[i][k]*a[i][k];
					}
					f			=	a[i][l];
					g			= -sign(sqrt(h),f);
					e[i]		=	scale*g;
					h		  -=	f*g;
					a[i][l]	=	f-g;
					f			=	0.0;
					for (j=0;j<=l;j++) {
						a[j][i]	= a[i][j]/h;
						g			= 0.0;
						for (k=0;k<=j;k++)
							 g += a[j][k]*a[i][k];
						if (l>j)
							 for (k=j+1;k<=l;k++)
								 g += a[k][j]*a[i][k];
						e[j] = g/h;
						f	  += e[j]*a[i][j];
					}
					hh = f/(h+h);
					for (j=0;j<=l;j++) {
						f = a[i][j];
						g = e[j]-hh*f;
						e[j] = g;
						for (k=0;k<=j;k++)
						  a[j][k] -= f*e[k] + g*a[i][k];
					}
				}
			}
			else
				e[i] = a[i][l];
			d[i] = h;
		}
	}
	d[0] = 0.0;
	e[0] = 0.0;
	for (i=0;i<n;i++) {
		l = i-1;
		if (d[i]!=0.0) {
			for (j=0;j<=l;j++) {
				g = 0.0;
				for (k=0;k<=l;k++)
					g += a[i][k]*a[k][j];
				for (k=0;k<=l;k++)
					a[k][j] -= g*a[k][i];
			}
		}
		d[i] = a[i][i];
		a[i][i] = 1.0;
		if (l>=0) {
			for (j=0;j<=l;j++)
				a[i][j] = a[j][i]  = 0.0;
		}
	}
}


void tqli (float *d, float *e, int n, float z[3][3]) {
	int m,l,iter,i,k,ok;
	float s,r,p,g,f,dd,c,b;
	if (n>1) {
		for (i=1;i<n;i++)
			e[i-1] = e[i];
		e[n-1] = 0;
		for (l=0;l<n;l++) {
			iter = 0;
			do {
				m = l;
				ok = 1;
				while (ok && m<n-1)
				{
					dd = fabs(d[m]) + fabs(d[m+1]);
					ok = (dd!=dd+fabs(e[m]));
					if (ok)
						m++;
				}
				if (m!=l) {
					g = (d[l+1]-d[l])/(2.0*e[l]);
					r = sqrt(1.0+g*g);
					g = d[m]-d[l]+e[l]/(g+sign(r,g));
					s = c = 1.0;
					p = 0.0;
					for (i=m-1;i>=l;i--) {
						f = s*e[i];
						b = c*e[i];
						if (fabs(f) >= fabs(g)) {
							c = g/f;
							r = sqrt (1 + c*c);
							e[i+1] = f*r;
							s = 1.0/r;
							c *= s;
						}
						else {
							s = f/g;
							r = sqrt (1 + s*s);
							e[i+1] = g*r;
							c = 1/r;
							s *= c;
						}
						g = d[i+1]-p;
						r = (d[i]-g)*s + 2.0*c*b;
						p = s*r;
						d[i+1] = g+p;
						g = c*r-b;
						for (k=0;k<n;k++) {
							f = z[k][i+1];
							z[k][i+1] = s*z[k][i] + c*f;
							z[k][i  ] = c*z[k][i] - s*f;
						}
					}
					d[l] -= p;
					e[l]	= g;
					e[m]	= 0.0;
				}
			} while (m!=l && ++iter<=30);
			if (m!=l) {
				fprintf (stderr, "ERROR (tqli:strnsub.c): More than 30 iterations.\n");
				exit(1);
			}
		}
	}
}


void pure (float *qxx, float *qxy, float *qxz,
			  float *qyy, float *qyz, float *qzz) {
	float a[3][3], b[3][3];
	float d[3],e[3];
	int	i,j;
	/*  CALCULATE EIGENVECTORS  */
	a[0][0] = *qxx;
	a[0][1] = *qxy;
	a[0][2] = *qxz;
	a[1][0] = *qxy;
	a[1][1] = *qyy;
	a[1][2] = *qyz;
	a[2][0] = *qxz;
	a[2][1] = *qyz;
	a[2][2] = *qzz;
	tred2 (a, 3, d, e);
	tqli	(d, e, 3, a);
	/*  FORM STRAIN MATRIX	*/
	for (i=0;i<3;i++)
		if (d[i]<0.0)
			d[i] = -sqrt(-d[i]);
		else
			d[i] =  sqrt( d[i]);
	for (i=0;i<3;i++)
	for (j=0;j<3;j++)
		b[i][j] =	  a[i][0] * d[0] * a[j][0]
						+ a[i][1] * d[1] * a[j][1]
						+ a[i][2] * d[2] * a[j][2];
	*qxx = b[0][0];
	*qxy = b[0][1];
	*qxz = b[0][2];
	*qyy = b[1][1];
	*qyz = b[1][2];
	*qzz = b[2][2];
}

void eigenstrain (float qxx, float qxy, float qxz,
						float qyy, float qyz, float qzz,
						float *v1, float *v2, float *v3) {

	float a[3][3];
	float d[3], e[3];
	int	i;
	/*  CALCULATE EIGENVECTORS  */
	a[1][1] = qxx;
	a[1][2] = qxy;
	a[1][3] = qxz;
	a[2][1] = qxy;
	a[2][2] = qyy;
	a[2][3] = qyz;
	a[3][1] = qxz;
	a[3][2] = qyz;
	a[3][3] = qzz;
	tred2 (a, 3, d, e);
	tqli	(d, e, 3, a);
	/*  REDUCE EIGENVALUES	*/
	for (i=0;i<3;i++)
		if (d[i]<0.0)
			d[i] = -sqrt(-d[i]);
		else
			d[i] =  sqrt( d[i]);
	/*  PASS RESULTS TO EIGENVECTOR ARRAYS  */
	v1[0] = d[0];
	v2[0] = d[1];
	v3[0] = d[2];
	for (i=1;i<4;i++) {
		v1[i] = a[0][i-1];
		v2[i] = a[1][i-1];
		v3[i] = a[2][i-1];
	}
}

/*

void main () {
	float a[3][3];
	float qxx,qxy,qxz,qyy,qyz,qzz;
	qxx = qyy = qzz = 6.0;
	qxy = qyz = qxz = 5.0;

	a[0][0] = qxx;
	a[1][1] = qyy;
	a[2][2] = qzz;
	a[0][1] = a[1][0] = qxy;
	a[1][2] = a[2][1] = qyz;
	a[0][2] = a[2][0] = qxz;
	printf ("\nOriginal Matrix\n");
	writematrix (a, 3,3 );

	pure (&qxx,&qxy,&qxz, &qyy,&qyz, &qzz);

	a[0][0] = qxx;
	a[1][1] = qyy;
	a[2][2] = qzz;
	a[0][1] = a[1][0] = qxy;
	a[1][2] = a[2][1] = qyz;
	a[0][2] = a[2][0] = qxz;
	printf ("\n\"Pure\" Matrix\n");
	writematrix (a, 3,3 );

	inverse (&qxx,&qxy,&qxz, &qyy,&qyz, &qzz);

	a[0][0] = qxx;
	a[1][1] = qyy;
	a[2][2] = qzz;
	a[0][1] = a[1][0] = qxy;
	a[1][2] = a[2][1] = qyz;
	a[0][2] = a[2][0] = qxz;
	printf ("\nInverse Matrix\n");
	writematrix (a, 3,3 );

}

*/