cditemp.c

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

/*
************************************************************************
Include Files
************************************************************************
*/
#include <stdlib.h>
#include <stdio.h>
#include <math.h>
#include "cdhouse.h"
#include "cdalloc.h"
#include "cditemp.h"
#include "parse.h"

/*
************************************************************************
Defines
************************************************************************
*/
#define X 0
#define Y 1
#define Z 2
#define NDIR 3

#define BOOLEAN int
#define FALSE	 0
#define TRUE	 1


/*
************************************************************************
Type Definitions
************************************************************************
*/

/*
************************************************************************
Local Function Prototypes
************************************************************************
*/
void CalcLinearVelocity
	(
	double (*Velocity)[NDIR],
	double  *Mass,
	SEL_T   *Select,
	int	  NumPart,
	double  CenterOfMassVelocity[NDIR]
	);

void CalcAngularVelocity
	(
	double (*Coord)[NDIR],
	double (*Velocity)[NDIR],
	double  *Mass,
	SEL_T  *Select,
	int	  NumPart,
	double  AngVelocity [NDIR],
	double  CenterOfMass[NDIR]
	);

void AddVelocity
	(
	double (*Velocity)[NDIR],
	SEL_T	  *Select,
	int		NumPart,
	double	NewVelocity [NDIR]
	);

void AddAngularVelocity
	(
	double (*Coord   )[NDIR],
	double (*Velocity)[NDIR],
	SEL_T	  *Select,
	int		NumPart,
	double	AngVelocity [NDIR],
	double	CenterOfMass[NDIR]
	);


/*
************************************************************************
Exported Function
************************************************************************
*/

/*
Initialize temperature for all particles and in specified directions
*/
void InitialTemp
(
Particle_t *a,
Simulation_t *s,
double InputTemp,
BOOLEAN UseSelect,
BOOLEAN UseOldVelocity,
int xflag,
int yflag,
int zflag
)
	{
	int	 i;
	int	 ipart;
	int    index;
	int	 idir;
	int	 ndir;
	int	 nfree;
	double (*Velocity)[NDIR] = NULL;
	double (*Coord   )[NDIR] = NULL;
	double  *Mass = NULL;
	BOOLEAN DirectionFlag[NDIR];
	double  EnergySum;
	double  CurrentKinetic;
	double  TargetKinetic;
	double  VelocityFactor;
	double  CurTemp;
	double  ZeroVelocity[NDIR] = {0.0, 0.0, 0.0};
	int     TempIndex;

	/*  Test for Zero Mass	and  Number of Free Particles  */
	nfree = 0;
	LOOP(ipart, a->np)
	if 
	(
	!IS_FIX(ipart) && 
	(!UseSelect || IS_SELECT(ipart))
	)
		{
		if (a->mass[ipart]==0.0)
			{
			printf ("ERROR INITIALIZING TEMPERATURE.\n");
			printf ("    Particle %i has zero mass.\n", ipart);
			CleanAfterError();
			}
		nfree++;
		}

	/*  Print error message if no free particles  */
	if (nfree==0)
		{
		if (UseSelect)
			printf 
			("ERROR (itemp):  No free particles selected (all are fixed).\n");
		else
			printf ("ERROR (itemp):  No free particles (all are fixed).\n");
		CleanAfterError();
		}

	/*  Insure velocity is allocated  */
	if (a->v==NULL)
		ALLOCATE (a->v, double, NDIR*a->NumPartAlloc)

	/*  Test for temperature equal to zero  */
	if (InputTemp==0)
		{
		LOOP (i,a->np)
		if (!UseSelect || IS_SELECT(i))
			{
			a->v[i*NDIR+X] = 0.0;
			a->v[i*NDIR+Y] = 0.0;
			a->v[i*NDIR+Z] = 0.0;
			}
		return;
		}

	/*  Set Direction Flag	*/
	DirectionFlag[X] = xflag;
	DirectionFlag[Y] = yflag;
	DirectionFlag[Z] = zflag;

	/*  Sum Number of Directions	*/
	ndir = 0;
	LOOP (idir,NDIR)
		if (DirectionFlag[idir])
			ndir++;

	if (ndir==0)
		{
		printf ("ERROR:  No directions included for random velocities.\n");
		CleanAfterError();
		}


	/* Allocate room for working version of velocity and masses  */
	ALLOCATE (Coord,    Coord_t, nfree)
	ALLOCATE (Velocity, Coord_t, nfree)
	ALLOCATE (Mass,     double,  nfree)

	/*
	Copy velocities and masses from particle structure
	to temporary storage
	*/
	TempIndex=0;
	LOOP (ipart, a->np)
		{
		if 
		(
		!IS_FIX(ipart) &&
		(!UseSelect || IS_SELECT(ipart))
		)
			{
			Coord 	[TempIndex] [X] = a->cur [NDIR*ipart + X];
			Coord 	[TempIndex] [Y] = a->cur [NDIR*ipart + Y];
			Coord 	[TempIndex] [Z] = a->cur [NDIR*ipart + Z];
			Mass		[TempIndex] 	 = a->mass[ 	 ipart	 ];
			TempIndex++;
			}
		}

	/*  Insure number of free particles hasn't changed  */
	ASSERT(TempIndex==nfree)

	/*  Set trial velocities  */
	LOOP (idir, NDIR)
		{

		/*  Assign random velocities	*/
		if (DirectionFlag[idir])
			{
			LOOP (i, nfree)
				{
				if (UseOldVelocity)
					Velocity[i][idir] = randfloat();
				else
					Velocity[i][idir] = nrand()/sqrt(Mass[i]);
				}
			}

		/*  Assign Zero Velocity  */
		else
			{
			LOOP (i, nfree)
				Velocity[i][idir] = 0.0;
			}
		}

	/*  Zero Linear Momentum  */
	SetVelocity(Velocity, Mass, NULL, nfree, ZeroVelocity);

	/*  Zero Angular Momentum	*/
	SetAngularVelocity 
		(Coord, Velocity, Mass, NULL, nfree, ZeroVelocity);

	/*  Re-initialize Velocities to Zero if not randomized  */
	LOOP (idir, NDIR)
		if (!DirectionFlag[idir])
			LOOP (ipart, nfree)
				Velocity[ipart][idir] = 0.0;

	/*  Scale Kinetic Energy to Temperature  */
	EnergySum = 0.0;
	LOOP (idir, NDIR)
	LOOP (ipart, nfree)
		EnergySum += Mass[ipart] * Velocity[ipart][idir] * Velocity[ipart][idir];

	/*  Test integrity of routine  */
	ASSERT(EnergySum!=0.0)


	/*  Normalize kinetic energy	*/
	CurrentKinetic  = 0.5 * EnergySum;
	TargetKinetic	 = 0.5 * ndir * nfree * BK * InputTemp;
	VelocityFactor  = s->dtime * sqrt (TargetKinetic / CurrentKinetic);
	LOOP (ipart, nfree)
	LOOP (idir,  NDIR )
		Velocity[ipart][idir] *= VelocityFactor;

	/*  Store kinetic energy and temperature	*/
	a->ekin = TargetKinetic;
	a->temp = InputTemp;

	/*  Store new velocities and higher derivatives  */
	TempIndex=0;
	LOOP (ipart, a->np)
		{
		if 
		(
		!IS_FIX(ipart) &&
		(!UseSelect || IS_SELECT(ipart))
		)
			{
			LOOP (idir, NDIR)
				{
				index = ipart * NDIR + idir;
				a->v[index] = Velocity[TempIndex][idir];
				/*  Reset higher derivative (for Gear integration)  */
				if (a->c2!=NULL)  a->c2[index] = 0.0;
				if (a->c3!=NULL)  a->c3[index] = 0.0;
				if (a->c4!=NULL)  a->c4[index] = 0.0;
				if (a->c5!=NULL)  a->c5[index] = 0.0;
				}
			TempIndex++;
			}
		}

	/*  Free allocation of temporary arrays  */
	FREE(Mass)
	FREE(Velocity)
	FREE(Coord)

	/*  Re-scale velocities to account for any constraint  */
	if (!a->cons) return;
	ApplyConstraint(a);
	CurTemp = GetTemperature(a,s,UseSelect);
	if (CurTemp==0)
		return;
	VelocityFactor = sqrt(InputTemp/CurTemp);
	LOOP (ipart, a->np)
		{
		if 
		(
		!IS_FIX(ipart) &&
		(!UseSelect || IS_SELECT(ipart))
		)
			{
			a->v[NDIR*ipart+X] *= VelocityFactor;
			a->v[NDIR*ipart+Y] *= VelocityFactor;
			a->v[NDIR*ipart+Z] *= VelocityFactor;
			}
		}

	CalcKineticEnergy(a,s);
	}



/*
For selected particles, reset center of mass velocity
*/
void SetVelocity
(
double (*Velocity)[NDIR],
double  *Mass,
SEL_T	  *Select,
int	  NumPart,
double  NewVelocity[NDIR]
)
	{
	double OldVelocity[NDIR];
	double VelocityChange[NDIR];

	/*   Calculate previous center of mass velocity  */
	CalcLinearVelocity
		(Velocity, Mass, Select, NumPart, OldVelocity);

	/*  Add new, remove old velocity  (convert velocity) */
	VelocityChange[X] = NewVelocity[X] - OldVelocity[X];
	VelocityChange[Y] = NewVelocity[Y] - OldVelocity[Y];
	VelocityChange[Z] = NewVelocity[Z] - OldVelocity[Z];
	AddVelocity (Velocity, Select, NumPart, VelocityChange);

	}