cdalloc.c

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

/*
************************************************************************
File Includes
************************************************************************
*/
#include <stdlib.h>
#include <stdio.h>
#include <string.h>
#include "cdhouse.h"
#include "cdalloc.h"


/*
************************************************************************
Defines
************************************************************************
*/
#define NEIGH_ALLOC_INC 16384
#define PART_ALLOC_INC 1024

/*
************************************************************************
Macros
************************************************************************
*/
/*  Macros used in CopyParticle()  */

#define COPY_COORD(ARRAY,TARGET,SOURCE) \
	if (ARRAY != NULL) \
		memcpy  \
		(&ARRAY[NDIR*TARGET], &ARRAY[NDIR*SOURCE], NDIR*sizeof(ARRAY[0]));


#define COPY_ELEM(ARRAY,TARGET,SOURCE) \
	if (ARRAY != NULL) \
		memcpy (&ARRAY[TARGET], &ARRAY[SOURCE], sizeof(ARRAY[0]));




/*
************************************************************************
Local Function Prototypes
************************************************************************
*/
void ReAllocateExternalVectorList
(
ExternalVector_t ***List,
int NumPart,
int OldNumPart
);

/*
************************************************************************
Exported Functions
************************************************************************
*/


/*
Allocate particle required arrays.
Set optional array pointers to NULL
*/
Particle_t *palloc (int NumPartAlloc, int NumNeighAlloc)
	{
	int i;
	int iclamp;
	Particle_t *a = NULL;


	/*  Allocate particle structure  */
	ALLOCATE (a,							  Particle_t,		  1			  )

	/*  Set allocation value  */
	a->NumPartAlloc = NumPartAlloc;
	a->NumNeighAlloc = NumNeighAlloc;

	/*  Allocate required coordinate arrays  */
	if (NumPartAlloc!=0)
		{
		ALLOCATE (a->cur,                  double,      NDIR*NumPartAlloc)
	
		/*  Allocate other particle arrays  */
		ALLOCATE (a->NeighborListIndex,	  int,				  NumPartAlloc)
		ALLOCATE (a->NeighborListLength,   int,				  NumPartAlloc)
		ALLOCATE (a->Selection,            SEL_T,            NumPartAlloc)
		ALLOCATE (a->type,					  BYTE,				  NumPartAlloc)
		ALLOCATE (a->mass,					  double,			  NumPartAlloc)
		ALLOCATE (a->eatom,					  double,			  NumPartAlloc)
		}

	/*  Allocate neighbor list  */
	if (NumNeighAlloc!=0)
		ALLOCATE (a->NeighborList, 		  int,   			  NumNeighAlloc)

	/* Set optional array pointers to NULL  */
	a->ref	 = NULL;
	a->ngh	 = NULL;
	a->v		 = NULL;
	a->c2 	 = NULL;
	a->c3 	 = NULL;
	a->c4 	 = NULL;
	a->c5 	 = NULL;
	a->rtype  = NULL;
	a->cons	 = NULL;
	a->disp	 = NULL;
	a->Damp   = NULL;

	/*
	Set counters
	*/
	a->nfix	 = 0;

	/* INITIALIZE FLAGS */
	a->IsInitializedBox       = FALSE;
	a->IsInitializedMass      = FALSE;
	a->IsInitializedCoord     = FALSE;
	a->IsInitializedPotential = FALSE;
	a->IsInitializedSurface   = FALSE;

#if 0
	a->coordflag = FALSE;
	a->boxflag	 = FALSE;
#endif

	a->InvalidNeighborList = TRUE;
	LOOP (i, NDIR)
		a->surf[i] = TRUE;

	/*  Initialize counts and averages	*/
	a->ebath   =  0.0;
	a->np 	  =  0;
	a->ng 	  =  0;
	a->nsel	  =  0;
	a->nbin	  = 100;
	a->step	  =  0;
	a->run	  =  0;

	/*  Intialize title to default string *title*  */
	LOOP (i, 8)
		strcpy (a->title[i], "*title*");

	/*  Initialize Damping   */
	a->UseDamp = FALSE;

	/*  Initialize Cavity  */
	a->UseCavity = FALSE;

	/*  Initialize Pressure settings  */
	a->BoxMotionAlgorithm = BMA_NONE;
	a->BoxMotionGeometry  = BMG_ISOTROPIC;
   a->Pressure[X]        = 0.0;
   a->Pressure[Y]        = 0.0;
   a->Pressure[Z]        = 0.0;

	/*  Volume Clamp Info  */
	a->VolClampStep = 33.0;

	/*  Initialize Repeating Boundary Geometry  */
	a->BoundaryType = BD_ORTHORHOMBIC;

	/*  Clamp info  */
	LOOP (iclamp, NUM_CLAMP_TAG)
		{
		a->UseClamp [iclamp] = FALSE;
		a->ClampTemp[iclamp] = 0.0;
		a->ClampStep[iclamp] = 33.0;
		}



	/*  RETURN POINTER  */
	return (a);
	}


/*  FREE ALLOCATION FOR PARTICLE STRUCTURE AND ALL ITS ARRAYS	*/
void pfree (Particle_t *a)
	{
	if (a)
		{
		/*  FREE COORDINATE TYPE ARRAYS	*/
		FREE(a->cur )
		FREE(a->ngh )
		FREE(a->ref )
		FREE(a->v   )
		FREE(a->f   )
		FREE(a->c2  )
		FREE(a->c3  )
		FREE(a->c4  )
		FREE(a->c5  )
		FREE(a->disp)
		FREE(a->Damp)

		/*  FREE NORMAL ARRAYS	*/
		FREE(a->type)
		FREE(a->rtype)
		FREE(a->mass)
		FREE(a->eatom)
		FREE(a->NeighborListIndex)
		FREE(a->NeighborListLength)
		FREE(a->NeighborList)
		FREE(a->Selection)
		FREE(a->tag)

		/*  Set counters to zero  */
		a->nfix	 = 0;

		/*  NOTE:  At some future time add deallocation of &a->cons  */
		}
	}


/*
Re-allocate or Newly Allocate particle coordinates and
all arrays whose size depends on number of particles.
(if they are already in use, determined if pointer not NULL)
Arrays such as neighbor list arrays are not affected.
*/
void ReallocateParticle (Particle_t *a, int NewNumPart)
	{
	int OldNumCoord;
	int OldNumPart;
	int NewAddPartAlloc;
	int NewAllocCoord;
	int NewAllocPart;
	int ipart;

	/*  Return if particles already allocated  */
	if (NewNumPart <= a->NumPartAlloc)
		return;

	/*  Determine number of additional particle allocations */
	NewAddPartAlloc = NewNumPart - a->NumPartAlloc;
	
	/*  Don't allocate less than PART_ALLOC_INC  */
	if (NewAddPartAlloc>0 && NewAddPartAlloc<PART_ALLOC_INC)
		{
		NewAddPartAlloc = PART_ALLOC_INC;
		}
	NewAllocPart = a->NumPartAlloc + NewAddPartAlloc;

	/*  Determine number of coordinates (3 times particles)  */
	OldNumPart = a->np;
	OldNumCoord   =  NDIR * OldNumPart;
	NewAllocCoord =  NDIR * NewAllocPart;


	/*  Allocate required coordinate arrays  */
	REALLOC(a->cur,                double, OldNumCoord, NewAllocCoord)

	/*  Allocate other required arrays  */
	REALLOC(a->NeighborListIndex,  int,    OldNumPart , NewAllocPart )
	REALLOC(a->NeighborListLength, int,    OldNumPart , NewAllocPart )
	REALLOC(a->Selection,          SEL_T,  OldNumPart , NewAllocPart )
	REALLOC(a->tag,                BYTE,   OldNumPart , NewAllocPart )
	REALLOC(a->type,               BYTE,   OldNumPart , NewAllocPart )
	REALLOC(a->mass,               double, OldNumPart , NewAllocPart )
	REALLOC(a->eatom,              double, OldNumPart , NewAllocPart )

	/*  Reallocate optional coordinate arrays  */
	if (a->ngh   != NULL)   REALLOC(a->ngh,  double,OldNumCoord,NewAllocCoord)
	if (a->ref   != NULL)   REALLOC(a->ref,  double,OldNumCoord,NewAllocCoord)
	if (a->v     != NULL)   REALLOC(a->v  ,  double,OldNumCoord,NewAllocCoord)
	if (a->f     != NULL)   REALLOC(a->f  ,  double,OldNumCoord,NewAllocCoord)
	if (a->c2    != NULL)   REALLOC(a->c2 ,  double,OldNumCoord,NewAllocCoord)
	if (a->c3    != NULL)   REALLOC(a->c3 ,  double,OldNumCoord,NewAllocCoord)
	if (a->c4    != NULL)   REALLOC(a->c4 ,  double,OldNumCoord,NewAllocCoord)
	if (a->c5    != NULL)   REALLOC(a->c5 ,  double,OldNumCoord,NewAllocCoord)
	if (a->disp  != NULL)   REALLOC(a->disp, double,OldNumCoord,NewAllocCoord)
	if (a->Damp  != NULL)   REALLOC(a->Damp, double,OldNumCoord,NewAllocCoord)

	/*  Reallocate other arrays  */
	if (a->rtype != NULL)   REALLOC(a->rtype,BYTE,OldNumPart,NewAllocPart)

	/*  For constraint array, free unneeded nodes if array will shrink  */
	if (NewAllocPart < OldNumPart && a->cons != NULL)
		for (ipart=NewNumPart; ipart<OldNumPart; ipart++)
			{
			FREE (a->cons[ipart])
			}

	/*  Now reallocate cons array  */
	if (a->cons  != NULL)	
		REALLOC (a->cons,Constraint_t*, OldNumPart, NewAllocPart)

	/*
	Reallocate External Vector arrays (SpringPtrList and ForcePtrList)
	*/
	ReAllocateExternalVectorList (&a->SpringPtrList, NewAllocPart, OldNumPart);
	ReAllocateExternalVectorList (&a->ForcePtrList,  NewAllocPart, OldNumPart);


	/*  Save new number of particles  */
#if 0
/* 26 March 1998  */
	a->np = NewNumPart;
#endif
	a->NumPartAlloc = NewAllocPart;

	}


/*  
New neighbor allocation routine
*/
void ReallocateNeighbors(Particle_t *a, int NewNumNeighbors)
	{
	WORD NewAllocation;
	WORD OldAllocation;
	WORD NewAllocationRounded;
	
	/*  No new neighbors needed  */
	if (NewNumNeighbors <= a->NumNeighAlloc)
		return;

	/*  Calculate new allocation  */
	OldAllocation = a->NumNeighAlloc;
	NewAllocation = a->NumNeighAlloc + NewNumNeighbors;

	/*  
	Round up New Neighbor allocation in units of NEIGH_ALLOC_INC  
	*/
	NewAllocationRounded = 
		NEIGH_ALLOC_INC * 
		( (NewAllocation-1)/NEIGH_ALLOC_INC + 1);

	ASSERT(NewAllocationRounded>=NewAllocation)

	/*  Reallocation neighbor list memory  */
	REALLOC(a->NeighborList,  int,    OldAllocation, NewAllocationRounded)

	/*  Record new allocation  */
	a->NumNeighAlloc = NewAllocationRounded;
	}



/*
Copy all relevant information from one particle to another
*/
void CopyParticle (Particle_t *a, int Target, int Source)
	{
#ifdef OLD_COPY
	int TargC;
	int SourC;
	int SizeC;

	TargC = NDIR * Target;
	SourC = NDIR * Source;
	SizeC = NDIR * sizeof (a->cur[0]);
#endif

	/*  Copy coordinate type information  */
	COPY_COORD (a->cur , Target,  Source)
	COPY_COORD (a->ngh , Target,  Source)
	COPY_COORD (a->v   , Target,  Source)
	COPY_COORD (a->ref , Target,  Source)
	COPY_COORD (a->f   , Target,  Source)
	COPY_COORD (a->c2  , Target,  Source)
	COPY_COORD (a->c3  , Target,  Source)
	COPY_COORD (a->c4  , Target,  Source)
	COPY_COORD (a->c5  , Target,  Source)
	COPY_COORD (a->disp, Target,  Source)
	COPY_COORD (a->Damp, Target,  Source)


#ifdef OLD_SEL
	COPY_ELEM  (a->sel,       Target, Source)
	COPY_ELEM  (a->fix,       Target, Source)
	COPY_ELEM  (a->set,       Target, Source)
#endif
	COPY_ELEM  (a->Selection, Target, Source)
	COPY_ELEM  (a->tag,       Target, Source)
	COPY_ELEM  (a->type ,     Target, Source)
	COPY_ELEM  (a->mass ,     Target, Source)
	COPY_ELEM  (a->eatom,     Target, Source)
	COPY_ELEM  (a->rtype,     Target, Source)


	/*  Copy constraint pointers (if array exists)	*/
	if (a->cons  != NULL)
		{

		/*  Free existing info at target  */
		FREE (a->cons[Target])

		/*  Copy Source pointer to Target  */
		a->cons[Target] = a->cons[Source];

		/*  Set value of Source pointer to null  */
		a->cons[Source] = NULL;
		}


	/*  Copy Spring list pointers (if array exists)  */
	if (a->SpringPtrList  != NULL)
		{

		/*  Free existing info at target  */
		FREE (a->SpringPtrList[Target])

		/*  Copy Source pointer to Target  */
		a->SpringPtrList[Target] = a->SpringPtrList[Source];

		/*  Set value of Source pointer to null  */
		a->SpringPtrList[Source] = NULL;
		}

	/*  Copy Spring list pointers (if array exists)  */
	if (a->ForcePtrList	!= NULL)
		{

		/*  Free existing info at target  */
		FREE(a->ForcePtrList[Target])

		/*  Copy Source pointer to Target  */
		a->ForcePtrList[Target] = a->ForcePtrList[Source];

		/*  Set value of Source pointer to null  */
		a->ForcePtrList[Source] = NULL;
		}

	}


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

void ReAllocateExternalVectorList
(
ExternalVector_t ***List,
int NumPart,
int OldNumPart
)
	{
	int ipart;

	/*  Return if no list  */
	if (*List==NULL)
		return;

	/*  Release storage for unneeded nodes  */
	if (NumPart<OldNumPart)
		for (ipart=NumPart; ipart<OldNumPart; ipart++)
			FREE((*List)[ipart])

	/*  Now reallocate pointer list	*/
	REALLOC (*List, ExternalVector_t*, OldNumPart, NumPart)

	}