direct.c

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

/*
************************************************************************
Compile Switches
************************************************************************
*/


/*
#define TEST_DRIVER
#define DEBUG
*/





/*
************************************************************************
Include Files
************************************************************************
*/
#include <stdio.h>
#include <stdlib.h>
#include <math.h>


/*
************************************************************************
Defines
************************************************************************
*/

#define NUM_FACTOR 3
#define X 0
#define Y 1
#define Z 2
#define NDIR 3
#define BOOLEAN int



/*
************************************************************************
Macros
************************************************************************
*/
#define LOOP(INDEX,LIMIT) \
	for ((INDEX)=0; (INDEX)<(LIMIT); (INDEX)++)


/*
Debugging Macro:  Print out variable name and value
*/
#define WRITEVAR(VAR_NAME,VAR_TYPE) \
      { \
      printf ("FILE %s LINE %i :", __FILE__, __LINE__); \
      printf ("%s = ", #VAR_NAME); \
      printf (#VAR_TYPE, (VAR_NAME) ); \
      printf ("\n"); \
      }

/*
Debugging Macro:  Print out debugging note
*/
#define WRITENOT(MSG) \
      { \
      printf ("FILE %s LINE %i :", __FILE__, __LINE__); \
      printf ("  (%s)\n", MSG); \
      }



/*
************************************************************************
Local Function Prototypes
************************************************************************
*/

int  GetCommonFactor (int[NDIR]);



/*
************************************************************************
Exported Function
************************************************************************
*/
void ScaleDirectionToMiller(double Dir[NDIR])
	{
	int i;
	int    IntDir  [NDIR];
	double SignDir [NDIR];
	double ScoreDir[NDIR];
	double Score;
	double Factor;
	double BestFitScore;
	double BestFitFactor;
	double NormalFactor;
	double CommonFactor;

	/*  Store direction signs, convert to abs()  */
	LOOP (i, NDIR)
		{
		if (Dir[i] < 0)
			SignDir[i] = -1;
		else
			SignDir[i] =  1;
		Dir[i] = fabs(Dir[i]);
		}

	/*  Calculate normal factor  */
	NormalFactor = 
		sqrt
			(
			Dir[X]*Dir[X] +
			Dir[Y]*Dir[Y] +
			Dir[Z]*Dir[Z] 
			);

	/*  Test for zero direction  */
	if (NormalFactor == 0)
		return;

	/*  Normalize working direction  */
	Dir[X] /= NormalFactor;
	Dir[Y] /= NormalFactor;
	Dir[Z] /= NormalFactor;

	/*  Initialize BestFitScore to worst possible score  */
	BestFitScore  = 3.0;
	BestFitFactor = 1.0;

	/*  Try all possible integer magnitude squared  */
	for (i=36; i>0;  i--)
		{

		/*  Factor that scales vector to magnitude squared  */
		Factor = sqrt(i);

		/*  
		Find fitting score
		- the distance from scaled direction to pure integer direction
		*/
		ScoreDir[X] = fmod(Factor*Dir[X], 1.0);
		ScoreDir[Y] = fmod(Factor*Dir[Y], 1.0);
		ScoreDir[Z] = fmod(Factor*Dir[Z], 1.0);

		/*  
		If difference between direction and integer is
		greater than 0.5, then reverse scale
		(for example, an error of 0.9 is really an error of 0.1)
		*/
		if (ScoreDir[X]>0.5) ScoreDir[X] = 1 - ScoreDir[X];
		if (ScoreDir[Y]>0.5) ScoreDir[Y] = 1 - ScoreDir[Y];
		if (ScoreDir[Z]>0.5) ScoreDir[Z] = 1 - ScoreDir[Z];

		/*
		Get final score
		*/
		Score = 
			ScoreDir[X]*ScoreDir[X] +
			ScoreDir[Y]*ScoreDir[Y] +
			ScoreDir[Z]*ScoreDir[Z] ;

		/*  
		Weight score (goodness of fit) so that larger
		vectors are discouraged
		*/
#if 0
		Score *= Factor;
#endif

#ifdef DEBUG
		/*  Print current factor and score  */
		printf ("i Factor Score: %5i %8.2lf %8.2lf\n", 
			i, Factor, Score);
#endif
#ifdef DEBUG
			/*  Print current best fit  */
		   printf ("i BestFactor BestScore: %5i %8.2lf %8.2lf\n", 
				i, BestFitFactor, BestFitScore);
#endif

		/*  Save current fitting info if best fit  */
		if (Score <= BestFitScore)
			{
			BestFitScore = Score;
			BestFitFactor= Factor;

			}
		}

	
	/*  Reconstruct absolute value direction vector  */
	Factor = BestFitFactor;
	Dir[X] *= Factor;
	Dir[Y] *= Factor;
	Dir[Z] *= Factor;

	/*  Factor out common integer factors  */
	IntDir[X] = Dir[X]+0.5;
	IntDir[Y] = Dir[Y]+0.5;
	IntDir[Z] = Dir[Z]+0.5;
	CommonFactor = GetCommonFactor (IntDir);

	/*  Reconstruct signed direction vector  */
	Dir[X] /= SignDir[X]*CommonFactor;
	Dir[Y] /= SignDir[Y]*CommonFactor;
	Dir[Z] /= SignDir[Z]*CommonFactor;

	}


/*
************************************************************************
Local Functions
************************************************************************
*/
/*
Remove common factors up to 5
*/
int GetCommonFactor (int InputDir[NDIR])
	{
	int FactorFound;
	int FactorTable[NUM_FACTOR] = {2, 3, 5};
	int TotalFactor = 1;
	int Dir[NDIR];
	int i;

	/*  Make local copy of directions  */
	Dir[X] = InputDir[X];
	Dir[Y] = InputDir[Y];
	Dir[Z] = InputDir[Z];

	/*  Test all factors from FactorTable[]  */
	LOOP (i, NUM_FACTOR)
		{

		/*  Divide out this factor repeatedly  */
		do
			{
			FactorFound = 
				(
				Dir[X]%FactorTable[i] == 0  &&
				Dir[Y]%FactorTable[i] == 0  &&
				Dir[Z]%FactorTable[i] == 0 
				);
			if (FactorFound)
				{
				Dir[X] /= FactorTable[i];
				Dir[Y] /= FactorTable[i];
				Dir[Z] /= FactorTable[i];
				TotalFactor *= FactorTable[i];
				}
			}
			while (FactorFound);
		}

	/*  Return largest common factor  */
	return TotalFactor;
	}



/*
************************************************************************
Test Driver
************************************************************************
*/


#ifdef TEST_DRIVER

int main (int argc, char *argv[])
	{
	double Dir[NDIR];

	if (argc<4)
		return 0;

	Dir[X] = atof(argv[Y]);
	Dir[Y] = atof(argv[Z]);
	Dir[Z] = atof(argv[NDIR]);
	
	ScaleDirectionToMiller(Dir);

	printf ("Best Fit Direction: %5.2f %5.2f %5.2f\n",
		Dir[X], Dir[Y], Dir[Z]);

	return 0;
	}

#endif