cditemp.c

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

/*
For selected particles, reset center of mass AngVelocity
*/
void SetAngularVelocity
(
double (*Coord)	[NDIR],
double (*Velocity)[NDIR],
double  *Mass,
SEL_T   *Select,
int	  NumPart,
double  NewAngVelocity[NDIR]
)
	{
	double OldAngVelocity[NDIR];
	double AngVelocityChange[NDIR];
	double CenterOfMass[NDIR];

	/*   Calculate previous center of mass AngVelocity  */
	CalcAngularVelocity
		(Coord, Velocity, Mass, Select, NumPart, OldAngVelocity, CenterOfMass);

	/*  Add new, remove old AngVelocity  */
	AngVelocityChange[X] = NewAngVelocity[X] - OldAngVelocity[X];
	AngVelocityChange[Y] = NewAngVelocity[Y] - OldAngVelocity[Y];
	AngVelocityChange[Z] = NewAngVelocity[Z] - OldAngVelocity[Z];
	AddAngularVelocity
		(Coord, Velocity, Select, NumPart, AngVelocityChange, CenterOfMass);

	}


/*
************************************************************************
Local Functions
************************************************************************
*/

void CalcLinearVelocity
(
double (*Velocity)[NDIR],
double *Mass,
SEL_T  *Select,
int	 NumPart,
double CenterOfMassVelocity[NDIR]
)
	{
	int ipart;
	int NumSelected;
	double TotalMass;

	/*  Initilize accumulators  */
	TotalMass = 0.0;
	CenterOfMassVelocity[X] = 0.0;
	CenterOfMassVelocity[Y] = 0.0;
	CenterOfMassVelocity[Z] = 0.0;


	/*  Sum velocity and mass	*/
	NumSelected = 0;
	LOOP(ipart,NumPart)
	if (Select==NULL || IS_SELECT2(Select,ipart))
		{
		CenterOfMassVelocity[X] 	+= Mass[ipart] * Velocity[ipart][X];
		CenterOfMassVelocity[Y] 	+= Mass[ipart] * Velocity[ipart][Y];
		CenterOfMassVelocity[Z] 	+= Mass[ipart] * Velocity[ipart][Z];
		TotalMass  += Mass[ipart];
		NumSelected++;
		}

	/*  If no atoms selected, then return zero velocity  */
	if (NumSelected==0)
		{
		CenterOfMassVelocity[X] = 0.0;
		CenterOfMassVelocity[Y] = 0.0;
		CenterOfMassVelocity[Z] = 0.0;
		return;
		}

	/*  Insure that total mass is not zero  */
	if (TotalMass==0.0)
		{
		printf ("ERROR IN TOTAL LINEAR MOMENTUM CALCULATION.\n");
		printf ("    Total mass is zero.\n");
		CleanAfterError();
		}

	/*  Normalize velocity	*/
	CenterOfMassVelocity[X] /= TotalMass;
	CenterOfMassVelocity[Y] /= TotalMass;
	CenterOfMassVelocity[Z] /= TotalMass;

	}



/*
Returns Angular Velocity about the Center of Mass,
and the Center of Mass
*/
void CalcAngularVelocity
(
double (*Coord)	  [NDIR],
double (*Velocity)  [NDIR],
double *Mass,
SEL_T  *Select,
int	 NumPart,
double AngVelocity  [NDIR],
double CenterOfMass [NDIR]
)
	{
	int	 ipart;
	int	 idir;
	int	 jdir;
	int	 i1,i2,i3, j1,j2,j3;
	int	 NumSelected;
	double AngMomentum[NDIR];
	double Lever	 [NDIR];
	double LeverSqr [NDIR];
	double MomentOfInertia	 [NDIR][NDIR];
	double InvMomentOfInertia[NDIR][NDIR];
	double TotalMass;
	double InvTotalMass;
	double Determinant;
	double MaxMatrixElement;
	double MinMatrixElement;


	/*   Calculate Center of Mass  */
	TotalMass = 0.0;
	CenterOfMass[X] = CenterOfMass[Y] = CenterOfMass[Z] = 0.0;
	NumSelected = 0;
	LOOP (ipart, NumPart)
	if (Select==NULL || IS_SELECT2(Select,ipart))
		{
		CenterOfMass[X] += Mass[ipart] * Coord[ipart][X];
		CenterOfMass[Y] += Mass[ipart] * Coord[ipart][Y];
		CenterOfMass[Z] += Mass[ipart] * Coord[ipart][Z];
		TotalMass  += Mass[ipart];
		NumSelected++;
		}

	/*  If no atoms selected, return zero angular velocity  */
	if (NumSelected==0)
		{
		AngVelocity[X] = 0.0;
		AngVelocity[Y] = 0.0;
		AngVelocity[Z] = 0.0;
		return;
		}

	/*  Check for zero mass  */
	if (TotalMass==0.0)
		{
		printf ("*** ERROR IN TOTAL ANGULAR MOMENTUM CALCULATION.\n");
		printf ("    Total mass is zero.\n");
		CleanAfterError();
		}

	/*  Normalize center of mass	*/
	InvTotalMass = 1.0/TotalMass;
	CenterOfMass[X] *= InvTotalMass;
	CenterOfMass[Y] *= InvTotalMass;
	CenterOfMass[Z] *= InvTotalMass;


	/*  Calculate Angular Momentum  */
	AngMomentum[X] = AngMomentum[Y] = AngMomentum[Z] = 0.0;
	LOOP (ipart, NumPart)
	if (Select==NULL || IS_SELECT2(Select,ipart))
		{
		Lever[X] = Coord[ipart][X] - CenterOfMass[X];
		Lever[Y] = Coord[ipart][Y] - CenterOfMass[Y];
		Lever[Z] = Coord[ipart][Z] - CenterOfMass[Z];
		AngMomentum[X] += Mass[ipart]*
			(Lever[Y]*Velocity[ipart][Z]-Lever[Z]*Velocity[ipart][Y]);
		AngMomentum[Y] += Mass[ipart]*
			(Lever[Z]*Velocity[ipart][X]-Lever[X]*Velocity[ipart][Z]);
		AngMomentum[Z] += Mass[ipart]*
			(Lever[X]*Velocity[ipart][Y]-Lever[Y]*Velocity[ipart][X]);
		}

	AngMomentum[X] *= InvTotalMass;
	AngMomentum[Y] *= InvTotalMass;
	AngMomentum[Z] *= InvTotalMass;


	/*  Calculate Moment of Inertia matrix  */
	LOOP (idir, NDIR)
	LOOP (jdir, NDIR)
		MomentOfInertia[idir][jdir] = 0.0;

	LOOP (ipart, NumPart)
	if (Select==NULL || IS_SELECT2(Select,ipart))
		{
		Lever[X]  = Coord[ipart][X] - CenterOfMass[X];
		Lever[Y]  = Coord[ipart][Y] - CenterOfMass[Y];
		Lever[Z]  = Coord[ipart][Z] - CenterOfMass[Z];
		LeverSqr[X] = Lever[X]*Lever[X];
		LeverSqr[Y] = Lever[Y]*Lever[Y];
		LeverSqr[Z] = Lever[Z]*Lever[Z];
		MomentOfInertia[X][X] += Mass[ipart] * (LeverSqr[Y]+LeverSqr[Z]);
		MomentOfInertia[Y][Y] += Mass[ipart] * (LeverSqr[Z]+LeverSqr[X]);
		MomentOfInertia[Z][Z] += Mass[ipart] * (LeverSqr[X]+LeverSqr[Y]);
		MomentOfInertia[Y][Z] -= Mass[ipart] *  Lever[Y]*Lever[Z];
		MomentOfInertia[Z][X] -= Mass[ipart] *  Lever[Z]*Lever[X];
		MomentOfInertia[X][Y] -= Mass[ipart] *  Lever[X]*Lever[Y];
		}
	MomentOfInertia[Z][Y] = MomentOfInertia[Y][Z];
	MomentOfInertia[X][Z] = MomentOfInertia[Z][X];
	MomentOfInertia[Y][X] = MomentOfInertia[X][Y];

	/*  Set very small values to zero  */
	MaxMatrixElement = 0.0;
	LOOP(idir, NDIR)
	LOOP(jdir, NDIR)
		if (fabs(MomentOfInertia[idir][jdir])>MaxMatrixElement)
			MaxMatrixElement = fabs(MomentOfInertia[idir][jdir]);

	MinMatrixElement = 1.0e-16 * MaxMatrixElement;
	LOOP(idir,NDIR)
	LOOP(jdir,NDIR)
		if (fabs(MomentOfInertia[idir][jdir])<MinMatrixElement)
			MomentOfInertia[idir][jdir] = 0.0;

	/*  Normalize	*/
	LOOP (idir, NDIR)
	LOOP (jdir, NDIR)
		MomentOfInertia[idir][jdir] *= InvTotalMass;

	/*
	Calculate transpose of Inverse of Moment of Inertia
		- not normalized by determinant
	*/
	LOOP(i1,NDIR)
		{
		i2 = (i1+1) % 3;
		i3 = (i1+2) % 3;
		LOOP(j1,NDIR)
			{
			j2 = (j1+1) % 3;
			j3 = (j1+2) % 3;
			InvMomentOfInertia[i1][j1] =
				MomentOfInertia[i2][j2]*MomentOfInertia[i3][j3] -
				MomentOfInertia[i2][j3]*MomentOfInertia[i3][j2];
			}
		}

	/*  Calculate determinant	*/
	Determinant = 0.0;
	LOOP(idir,NDIR)
	LOOP(jdir,NDIR)
		Determinant +=
			MomentOfInertia[idir][jdir]*InvMomentOfInertia[idir][jdir];
	Determinant /= 3.0;

	if (Determinant==0.0)
		{
		printf ("*** WARNING:  MOMENT OF INERTIA IS SINGULAR\n");
		printf ("    (mass distribution is confined to a plane,\n");
		printf ("     line or point).\n");
		IncrementNumberOfWarnings();
		return;
		}

	if (Determinant<0.0)
		{
		printf ("ERROR: Moment of inertia has negative determinant\n");
		printf ("Do you have 2 or fewer particles?\n");
		printf ("Or, are do all your particles lie on a single line?\n");
		CleanAfterError();
		}


	/*  Calculate Angular Velocity  */
	LOOP(idir,NDIR)
		{
		AngVelocity[idir] =
			(
			InvMomentOfInertia[idir][X]*AngMomentum[X] +
			InvMomentOfInertia[idir][Y]*AngMomentum[Y] +
			InvMomentOfInertia[idir][Z]*AngMomentum[Z]
			)
			/Determinant;
		}
	}




/*  Add velocity to all selected particles  */
void AddVelocity
(
double (*Velocity)[NDIR],
SEL_T  *Select,
int NumPart,
double CenterOfMassVelocity[NDIR]
)
	{
	int ipart;

	/*   Add velocity to selected particles  */
	LOOP (ipart, NumPart)
	if (Select==NULL || IS_SELECT2(Select,ipart))
		{
		Velocity[ipart][X] += CenterOfMassVelocity[X];
		Velocity[ipart][Y] += CenterOfMassVelocity[Y];
		Velocity[ipart][Z] += CenterOfMassVelocity[Z];
		}
	}



/*  Add angular velocity to selected particles	*/
void AddAngularVelocity
(
double (*Coord)	[NDIR],
double (*Velocity)[NDIR],
SEL_T	  *Select,
int		NumPart,
double	AngVelocity[NDIR],
double	CenterOfMass[NDIR]
)
	{
	int ipart;
	double Lever[NDIR];

	/*  Add Angular Velocity  */
	LOOP (ipart, NumPart)
	if (Select==NULL || IS_SELECT2(Select,ipart))
		{

		/*  Find position of particle relative to center of mass  */
		Lever[X] = Coord[ipart][X] - CenterOfMass[X];
		Lever[Y] = Coord[ipart][Y] - CenterOfMass[Y];
		Lever[Z] = Coord[ipart][Z] - CenterOfMass[Z];

		/*  Add linear velocity due to angular velocity  */
		Velocity[ipart][X] += AngVelocity[Y]*Lever[Z] - AngVelocity[Z]*Lever[Y];
		Velocity[ipart][Y] += AngVelocity[Z]*Lever[X] - AngVelocity[X]*Lever[Z];
		Velocity[ipart][Z] += AngVelocity[X]*Lever[Y] - AngVelocity[Y]*Lever[X];
		}

	}