cddisp.c

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

/*
************************************************************************
File Includes
************************************************************************
*/
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include "cdhouse.h"
#include "cdsubs.h"
#include "cdalloc.h"
#include "iomngr.h"
#include "strsub.h"
#include "parse.h"

/*
************************************************************************
Defines
************************************************************************
*/
#define DEFAULT_DISP_MOVE_SCALE 1.0
#define ANG_TO_CM 1e-8
#define CM_TO_ANG 1e+8

/*
************************************************************************
External Variables
************************************************************************
*/
extern Particle_t *a;
extern LIST 	 *inlist;


/*
************************************************************************
Module-Wide variables
************************************************************************
*/
BOOLEAN UseSelect_m;

/*
************************************************************************
Local Function Prototypes
************************************************************************
*/
void AllocateDispArray (void);
void ReadDisp_Clear	 (char *InputString);
void ReadDisp_Move	 (char *InputString);
void ReadDisp_Read	 (char *InputString);
void ReadDisp_RefCalc (char *InputString);
void ReadDisp_Scale	 (char *InputString);
void ReadDisp_Set 	 (char *InputString);


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


/*  Read disp command  */
void read_disp (char *InputString)
	{
	char *LeadingToken;

	/*  Get leading token  */
	LeadingToken = strhed (&InputString);

	/*  Test for SEL option  */
	UseSelect_m = !strcmpi (LeadingToken, "SEL");

	/*  For for SEL option but no particle selected  */
	if (UseSelect_m)
		CheckForNoneSelected();

	/*  If SEL option, get next token  */
	if (UseSelect_m)
		LeadingToken = strhed (&InputString);

	/*  Clear option	*/
	if 	  (!strcmpi (LeadingToken, "CLEAR"  ))
		ReadDisp_Clear   (InputString);

	/*  Move option  */
	if 	  (!strcmpi (LeadingToken, "MOVE"  ))
		ReadDisp_Move	  (InputString);

	/*  Read option  */
	else if (!strcmpi (LeadingToken, "READ"   ))
		ReadDisp_Read	  (InputString);

	/*  Refcalc option  */
	else if (!strcmpi (LeadingToken, "REFCALC"))
		ReadDisp_RefCalc (InputString);

	/*  Scale option	*/
	else if (!strcmpi (LeadingToken, "SCALE"  ))
		ReadDisp_Scale   (InputString);

	/*  Set option  */
	else if (!strcmpi (LeadingToken, "SET"    ))
		ReadDisp_Set	  (InputString);

	/*  Unknown option  */
	else
		{
		printf ("ERROR (read_disp):  Unknown option (%s) to DISP command.\n",
			LeadingToken);
		CleanAfterError();
		}

	}


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

/*
Allocate storage for disp array
*/
void AllocateDispArray (void)
	{
	/*  Allocate space for Disp array  */
	FREE(a->disp)
	ALLOCATE(a->disp, double, NDIR*a->NumPartAlloc)
	}

/*
Clear selected displacements to zero
(If all displacements selected, release storage)
*/
#pragma argsused
void ReadDisp_Clear	 (char *InputString)
	{
	int ipart;
	int NumPart;
	double (*Displacement)[NDIR];

	/*  Test to see if Displacements exist  */
	if (a->disp==NULL)
		{
		return;
		}

	/*  If clearing all particles, then simply release disp array	*/
	if (!UseSelect_m || a->nsel==a->np)
		{
		FREE(a->disp)
		return;
		}

	/*  Initialize variables  */
	NumPart = a->np;
	Displacement = (double (*)[NDIR]) a->disp;

	/*  Zero selected displacements	*/
	for (ipart=0; ipart<NumPart; ipart++)
		{
		if (!IS_SELECT(ipart) )
			{
			Displacement[ipart][X] = 0.0;
			Displacement[ipart][Y] = 0.0;
			Displacement[ipart][Z] = 0.0;
			}
		}

	}


/*
Move particles by adding displacement to them
*/
void ReadDisp_Move	 (char *InputString)
	{
	int	 NumPart;
	int	 ipart;
	double (*Coord 	  )[NDIR];
	double (*Displacement)[NDIR];
	double DispScale;

	/*  Test to see if Displacements exist  */
	if (a->disp==NULL)
		{
		return;
		}

	/*  Test for scale option	*/
	if (!IsBlank (InputString))
		DispScale = dblstrf (&InputString);
	else
		DispScale = DEFAULT_DISP_MOVE_SCALE;


	/*  Intialize variables  */
	NumPart = a->np;
	Coord 		 = (double (*)[NDIR]) a->cur ;
	Displacement = (double (*)[NDIR]) a->disp;

	/*  Displace selected particles	*/
	for (ipart=0; ipart<NumPart; ipart++)
		{
		if (!UseSelect_m || IS_SELECT(ipart) )
			{
			Coord[ipart][X] += DispScale*Displacement[ipart][X];
			Coord[ipart][Y] += DispScale*Displacement[ipart][Y];
			Coord[ipart][Z] += DispScale*Displacement[ipart][Z];
			}
		}

	}

/*  Read Displacements	*/
void ReadDisp_Read	 (char *InputString)
	{
	char BufferString[NSTR];
	char *TempString;
	int  *IxReadToPart = NULL;
	int  NumPart;
	int  NumSel;
	int  nread;
	int  iread;
	int  ipart;
	int  idir;
	double (*Displacement)[NDIR];


	/*  Allocate disp array if not present  */
	if (a->disp==NULL)
		AllocateDispArray();

	/*  Initialize number of particles	*/
	NumPart = a->np;
	NumSel  = a->nsel;
	Displacement = (double (*)[NDIR]) a->disp;

	/*  Read number of input Displacements  */
	nread = intstrf (&InputString);

	/*  Insure number to read equals number selected  */
	if (UseSelect_m && nread!=NumSel)
		{
		printf ("ERROR (ReadDisp_Read):  ");
		printf ("Number of expected input displacements (%i) different from\n",
			nread);
		printf ("   number of selected atoms (%i).\n", NumSel);
		CleanAfterError();
		}

	/*  Insure number to read equals number of atoms  */
	else if (!UseSelect_m && nread!=NumPart)
		{
		printf ("ERROR (ReadDisp_Read):  ");
		printf ("Number of expected input displacements (%i) different from\n",
			nread);
		printf ("   number of atoms (%i).\n", NumPart);
		CleanAfterError();
		}

	/*  Allocate storage for particle indices  */
	ALLOCATE (IxReadToPart, int, nread)

	/*  Store particle indices for reading displacements	*/
	iread = 0;
	for (ipart=0; ipart<NumPart; ipart++)
		{
		if (!UseSelect_m || IS_SELECT(ipart) )
			{
			if (iread>=nread)
				{
				printf ("INTERNAL ERROR (ReadDisp_Read):  ");
				printf ("Number of selected atoms is wrong.\n");
				CleanAfterError();
				}
			IxReadToPart[iread] = ipart;
			iread++;
			}
		}

	/*  Read displacements	*/
	iread = 0;
	idir	= X;
	while (	 iread < nread
			 && NULL!=m_gets_list(BufferString,NSTR,inlist))
		{
		TempString = BufferString;
		while (!IsBlank (TempString))
			{
			ipart = IxReadToPart [iread];
			Displacement[ipart][idir] = ANG_TO_CM * dblstrf (&TempString);
			if (idir==Z)
				{
				idir = X;
				iread++;
				}
			else
				idir++;
			}
		}

	/*  Check that all expected displacements have been read in  */
	if (iread < nread)
		{
		printf ("ERROR (ReadDisp_Read):  Number of displacments read is %i.\n", iread);
		printf ("   Number expected is %i.\n", nread);
		CleanAfterError();
		}

	/*  Free allocation	*/
	FREE (IxReadToPart)

	}

/*  Calculate displacements from reference particles	*/
#pragma argsused
void ReadDisp_RefCalc (char *InputString)
	{
	int	 NumPart;
	int	 ipart;
	double (*Coord)[NDIR];
	double (*Reference)[NDIR];
	double (*Displacement)[NDIR];

	/*  Allocate disp array if not present  */
	if (a->disp==NULL)
		AllocateDispArray();

	/*  Intialize variables  */
	NumPart = a->np;
	Coord 		 = (double (*)[NDIR]) a->cur ;
	Reference	 = (double (*)[NDIR]) a->ref ;
	Displacement = (double (*)[NDIR]) a->disp;

	/*  Test for existence of reference particles  */
	if  (Reference==NULL)
		{
		printf ("ERROR (ReadDisp_RefCalc):  There are no reference particles.\n");
		CleanAfterError();
		}

	/*  Calculate displacements  */
	for (ipart=0; ipart<NumPart; ipart++)
		{
		if (!UseSelect_m || IS_SELECT(ipart) )
			{
			Displacement[ipart][X] = Coord[ipart][X] - Reference[ipart][X];
			Displacement[ipart][Y] = Coord[ipart][Y] - Reference[ipart][Y];
			Displacement[ipart][Z] = Coord[ipart][Z] - Reference[ipart][Z];
			}
		}
	}

void ReadDisp_Scale	 (char *InputString)
	{
	int	 NumPart;
	int	 ipart;
	double Scale;
	double (*Displacement)[NDIR];

	/*  Allocate disp array if not present  */
	if (a->disp==NULL)
		AllocateDispArray();

	/*  Intialize variables  */
	NumPart = a->np;
	Displacement = (double (*)[NDIR]) a->disp;

	/*  Read value of scale  */
	Scale = dblstrf (&InputString);

	/*  Apply scale to displacements  */
	for (ipart=0; ipart<NumPart; ipart++)
		{
		if (!UseSelect_m || IS_SELECT(ipart) )
			{
			Displacement[ipart][X] *= Scale;
			Displacement[ipart][Y] *= Scale;
			Displacement[ipart][Z] *= Scale;
			}
		}
	}

void ReadDisp_Set 	 (char *InputString)
	{
	int	 NumPart;
	int	 ipart;
	double (*Displacement)[NDIR];
	double InputDisplacement[NDIR];

	/*  If no displacement array, then don't scale  */
	if (a->disp==NULL)
		return;

	/*  Intialize variables  */
	NumPart = a->np;
	Displacement = (double (*)[NDIR]) a->disp;

	/*  Read value of displacement  */
	InputDisplacement[X] = dblstrf (&InputString);
	InputDisplacement[Y] = dblstrf (&InputString);
	InputDisplacement[Z] = dblstrf (&InputString);

	/*  Set Displacements  */
	for (ipart=0; ipart<NumPart; ipart++)
		{
		if (!UseSelect_m || IS_SELECT(ipart) )
			{
			Displacement[ipart][X] = InputDisplacement[X];
			Displacement[ipart][Y] = InputDisplacement[Y];
			Displacement[ipart][Z] = InputDisplacement[Z];
			}
		}
	}