cdeam.c

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

/*
************************************************************************
History
************************************************************************
*/
/*
11 Aug 1997
   (1)  in _free_eam() use FREE() macro instead of free().  The macro
        calls free() but then sets pointer value to NULL.  Previously
        the pointer was not being set to NULL, and subsequent calls
        due to a POTENTIAL SET EAM 2 command would try to free() the
        non-NULL pointer.
26 Jan 1998
    (1)  Changed all calls to free() into FREE() macros
 5 Jun 1998
    (1)  Store last value of potential in p[e->n-1]
    (2)  When r=rcut have __efn() return value stored above
14 Sep 1998
    (1)  Add warning messages when reading potential (using
         either standard for FORMULA option) without NTABLE,
         X0 or XN.
    (2)  Correct error in GetAtomSubset() which omitted the
         last atom when using odd number of aatoms.  This error 
         would have no effect when the last atom in array
         was also the last atom in the neighbor list, in which
         case it wouldn't be needed because it was handled
         as a neighbor of previous atoms.
         
*/




/*
************************************************************************
Compiler Switches
************************************************************************
*/
/*  Use PC Assembly code for fabs() and poly() routine  */
#define ASM_CODE
#undef  ASM_CODE

/*  Use special DEBUG file to store diagnostic output  */
#define DEBUG
#undef  DEBUG

/*  Use new thread code (May 12, 1998)  */
#define USE_THREAD_CODE

/*  Dependent compiler switches  */
#ifndef USE_THREAD_CODE
#define NOT_USE_THREAD_CODE
#endif

/*  Insure that USE_THREAD_CODE and NOT_USE_THREAD_CODE are not both defined  */
#ifdef  USE_THREAD_CODE
#ifdef  NOT_USE_THREAD_CODE
#undef  NOT_USE_THREAD_CODE
#endif
#endif


/*  Substitute __efn __dfn with explicit calls to __fn  */
#define USE_FN_MACRO
#undef  USE_FN_MACRO


/*  
If using INTEL, call FPU initialization at 
beginning of each thread function.  May 12, 1998
*/
#ifdef INTEL
#define INIT_INTEL_FPU asm("finit");
#else
#define INIT_INTEL_FPU
#endif



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


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

/*  STRUCTURE DEFINITIONS	*/
typedef struct
{
	double x0,xn,dx;
	int	 n;
	double *a;
} EAMdata_t;


/*  Type for passing info to threads and function */
typedef struct
	{
	int  First;
	int  Num;
	}
	ThreadParam_t;



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

/*  Constants defining potential types  */
#define POT_NONE	0
#define POT_PAIR	1
#define POT_DENS	2
#define POT_EMBED 3

/*  Size of input string buffer	*/
#define NUM_BUFFER 256

/*  Number of columns of potential data to a line	*/
#define MAX_COL 5


/*
************************************************************************
Macros
************************************************************************
*/
#ifdef USE_FN_MACRO
#define __pfn(TYPE,DERIV,RADIUS)  __fn(&(Pair_m[TYPE]),DERIV,RADIUS)
#define __dfn(TYPE,DERIV,RADIUS)  __fn(&(Dens_m[TYPE]),DERIV,RADIUS)
#endif



/*
************************************************************************
Local Function Prototypes
************************************************************************
*/
void	 __parse_pot (int *, int *, char *, char **);
void	 __write_eamtype_bin  (FILE *fout, EAMdata_t *e);
void	 __read_eamtype_bin	 (FILE *fout, EAMdata_t *e);

#ifndef USE_FN_MACRO
double __pfn (int, int, double);
double __dfn (int, int, 	  double);
#endif

double __efn (int, int, 	  double);
double __fn  (EAMdata_t *, int, double);
void	 __write_potstate_bin  (char *fname);
void	 __read_potstate_bin   (char *fname);
void	 __write_eamtype_bin   (FILE *fout, EAMdata_t *e);
void	 __read_eamtype_bin	  (FILE *fout, EAMdata_t *e);
void	 __write_potstate_text (char *fname);
void	 __read_potstate_text  (char *fname);
void	 __write_eamtype_text  (FILE *fout, EAMdata_t *e);
void	 __read_eamtype_text   (FILE *fout, EAMdata_t *e);
void	 __free_eam (void);

EAMdata_t *GetPotPtr 		(int, int);
void		  WritePotential	(char *);
int		  GetFirstType 	(int);
int		  GetSecondType	(int);

static void	 CalcCutoff (void);

static void GetAtomSubset  (int, int, int *, int *);
static void ThreadCalcRho  (void *);
static void ThreadCalcForce(void *);


#ifdef ASM_CODE
double __poly (double, int, double *);
#endif

/*
************************************************************************
External Variables
************************************************************************
*/
/*  EXTERNAL VARIABLES	*/
extern Simulation_t *s;
extern LIST 		*inlist;


/*
************************************************************************
Module Variables
************************************************************************
*/
/*  GLOBAL VARIABLES  */
int	 NumType_m=0;
static EAMdata_t *Pair_m    = NULL;
static EAMdata_t *Embed_m   = NULL;
static EAMdata_t *Dens_m    = NULL;
static double	*Rc_m 	    = NULL;
static double	*RcSqr_m     = NULL;
static double	*DensRc_m 	 = NULL;
static double	*DensRcSqr_m = NULL;
static double  *Rho_m       = NULL;
static double  *Fdrv_m      = NULL;
static Particle_t  *a_m     = NULL;
static double  Epair_m      = 0.0;

/*  Thread Variables  */
#ifdef USE_THREAD_CODE
static PTHREAD_MUTEX_T RhoLock_m;
static PTHREAD_MUTEX_T ForceLock_m;
static PTHREAD_MUTEX_T StressLock_m;
#endif



/*
************************************************************************
Macros
************************************************************************
*/

/*  Redefine fabs() function	*/
#ifdef ASM_CODE
double fabs_new(double);
#define FABS fabs_new
#else
#define FABS fabs
#endif


/*
************************************************************************
Exported Subroutines
************************************************************************
*/

#pragma argsused
void em_energy_list	(Particle_t *a, Simulation_t *s, int sflag)
	{
	unsigned char *type, t1, t2, tt;
	int i,j,k;
	double xd,yd,zd,xb,yb,zb,r,rsq,t;
	double xbox2,ybox2,zbox2, etemp, *eatom;
	double *c, *ck, *cj;
	double rc;
	double *rho = NULL;
	int xsurf = a->surf[X];
	int ysurf = a->surf[Y];
	int zsurf = a->surf[Z];
	int EndInterval;
	int StartInterval;
	EAMdata_t *e;


#ifdef DEBUG
	FILE *DebugFile;
	DebugFile = fopen ("Debug.fil", "wt");
	if (DebugFile==NULL)
		{
		printf ("ERROR:  Cannot open debug.fil\n");
		CleanAfterError();
		}
#endif

	/*  NEIGHBOR LIST FOR ALL PARTICLES  */
CHECK_HEAP
   a->CalcUniqueNeighbors = TRUE;
 
	search_all (a);
CHECK_HEAP

	/*  SET UP POINTER TO COORDINATES  */
	type	= a->type;
	eatom = a->eatom;

	/*  Allocate rho	*/
	ALLOCATE (rho, double, a->np);


	/*  INITIALIZE VALUES  */
	xb 	= a->bcur[X];
	yb 	= a->bcur[Y];
	zb 	= a->bcur[Z];
	xbox2 = xb - s->cutoff;
	ybox2 = yb - s->cutoff;
	zbox2 = zb - s->cutoff;

	/* INITIALIZE INDIVIDUAL ATOM ENERGIES  */
	/* CHECK TYPES 								 */
	for (i=0;i<a->np;i++)
		{
		a->eatom[i] = 0.0;
		if (type[i]>=NumType_m)
			{
			printf ("ERROR (em_energy_list):\n");
			printf ("  Type for particle %i\n", type[i]+1);
			printf ("  is out of range [1..%i]\n", NumType_m);
			CleanAfterError();
			}
		}


	/* INDIVIDUAL PAIR POTENTIALS  */
	c = a->cur;
	for (j=0; j<a->np; j++)
		{
		StartInterval = a->NeighborListIndex[j];
		EndInterval   = StartInterval + a->NeighborListLength[j];
		for (i=StartInterval; i<EndInterval; i++)
			{
			k	= a->NeighborList[i];
			cj = c + (j << 1) + j;
			ck = c + (k << 1) + k;
			t1 = type[j];
			t2 = type[k];
			tt = COMBINED_TYPE(t1,t2);
			rc = Rc_m[tt];
			/*  dislpacement relative to j  */
			xd = FABS(ck[0] - cj[0]);
			if (!xsurf)
				if (xd>xbox2)
					xd = xb - xd;
			if (xd<rc)
				{
				yd  = FABS(ck[1] - cj[1]);
				if (!ysurf)
					if (yd>ybox2)
						yd = yb - yd;
				if (yd<rc)
					{
					zd = FABS(ck[2] - cj[2]);
					if (!zsurf)
						if (zd>zbox2)
							zd = zb - zd;
					if (zd<rc)
						{
						rsq = xd*xd + yd*yd + zd*zd;
						/*  CALCULATE ENERGY CONTRIBUTIONS FOR CURRENT PAIR  */
						if (rsq<RcSqr_m[tt])
							{
							/*  CALCULATE ENERGY CONTRIBUTIONS FOR CURRENT PAIR  */
							r	= sqrt(rsq);
							if (t1==t2)
								{
								rho[j] += (t = __dfn(t2, 0, r));
								rho[k] +=  t;
								}
							else
								{
								rho[j]	+= __dfn (t2, 0, r);
								rho[k]	+= __dfn (t1, 0, r);
								}
							etemp 	 = __pfn (tt, 0, r);
#ifdef DEBUG
fprintf (DebugFile, "%i %i %14.4e %14.4e\n", t1, t2, r, etemp);
#endif
							eatom[j] += etemp;
							eatom[k] += etemp;
							}
						}
					}
				}
			}
		}

	/*  SUM POTENTIAL AND EMBEDDING ENERGIES	*/
	a->epot = 0.0;
	for (i=0;i<a->np;i++)
		{
		/*  Scale eatom by a half  */
		eatom[i] *= 0.5;

		/*  Calculate embedding function  */
		e = &Embed_m[type[i]];

		/*  No embedding function for this type  */
		if (e->a==NULL)
			{
			continue;
			}

		/*  Electron Density too high  */
		if (rho[i] > e->xn)
			{
			printf 
				("ERROR:  Electron density (%le) out of range [%le, %le]",
				rho[i], e->x0, e->xn);
			printf ("for atom %i (type %i).\n", i+1, type[i]+1);
			CleanAfterError();
			}

		/*  Calculate embedding and its derivative  */
		eatom[i] += __fn(e, 0, rho[i]);
		a->epot	+= eatom[i];
		}

	/*  FREE ARRAYS  */
	FREE (rho)

#ifdef DEBUG
	fclose (DebugFile);
#endif
	}


#pragma argsused
void em_read_potential (char *LeadingTokenPtr, char *InputString)
	{
	int PotType;
	int AtomType;
	int nread;
	int i,j,k,l,m,n;
	char InputBuffer[NUM_BUFFER];
	EAMdata_t *e = NULL;
	double *a = NULL;
	double *p = NULL;
	double *atmp = NULL;
	BOOLEAN ScaleInput;
	double  ScaleFactor;
	int	  itable;
	BOOLEAN FormulaOptionChosen;
	int	  NumFormula;
	int	  iformula;

	if (!strcmpi(LeadingTokenPtr,"INIT"))
		{
		int NumPair;

		/*  FREE POTENTIALS	*/
		__free_eam();

		/*  READ NEW POTENTIAL */
		NumType_m = intstrf (&InputString);
		if (NumType_m<=0)
			NumType_m = 0;

		/*  TEST FOR NumType_m==0	*/
		if (NumType_m==0)
			{
			printf ("ERROR (em_read_potential):  Must specify number of types.\n");
			CleanAfterError();
			}

		/*  ALLOCATE POTENTIALS  */
		NumPair	  = NumType_m*(NumType_m+1)/2;
		ALLOCATE (Rc_m,        double,    NumPair);
		ALLOCATE (RcSqr_m,     double,    NumPair);
		ALLOCATE (DensRc_m,    double,    NumPair);
		ALLOCATE (DensRcSqr_m, double,    NumPair);
		ALLOCATE (Dens_m,      EAMdata_t, NumType_m);
		ALLOCATE (Embed_m,     EAMdata_t, NumType_m);
		ALLOCATE (Pair_m,      EAMdata_t, NumPair);
		return;
		}

	else if (!strcmpi(LeadingTokenPtr,"WRITE_STATE_TEXT"))
		{
		__write_potstate_text (InputString);
		return;
		}

	else if (!strcmpi(LeadingTokenPtr,"READ_STATE_TEXT"))
		{
		__read_potstate_text  (InputString);
		return;
		}

	else if (!strcmpi(LeadingTokenPtr,"WRITE_STATE_BIN"))
		{
		__write_potstate_bin  (InputString);
		return;
		}

	else if (!strcmpi(LeadingTokenPtr,"READ_STATE_BIN"))
		{
		__read_potstate_bin	 (InputString);
		return;
		}

	else if (!strcmpi(LeadingTokenPtr,"EVAL"))
		{
		int nd;
		double x,f;

		LeadingTokenPtr = strhed (&InputString);
		__parse_pot (&PotType, &AtomType, LeadingTokenPtr, &InputString);
		nd = intstrf (&InputString);
		x	= dblstrf (&InputString);

		/* Evaluate potential and adjust units  */
		if (PotType==POT_PAIR)
			f = __pfn (AtomType, nd, 1e-8*x) * s->eunit;
		else if (PotType==POT_DENS)
			f = __dfn (AtomType, nd, 1e-8*x);
		else
			f = __efn (AtomType, nd,		x) * s->eunit;
		printf ("        %le\n", f);
		return;
		}

	/*  Scale potential	*/
	else if (!strcmpi(LeadingTokenPtr,"SCALE"))
		{
		/*  Get option - either "INPUT" or "OUTPUT"  */
		LeadingTokenPtr = strhed (&InputString);

		/*  Determine if scaling function input or output	*/
		if 	  (!strcmpi(LeadingTokenPtr, "INPUT") )
			ScaleInput = TRUE;
		else if (!strcmpi(LeadingTokenPtr, "OUTPUT"))
			ScaleInput = FALSE;
		else
			{
			printf
				("ERROR(em_read_potential):  Unrecognized option SCALE %s.\n",
				 LeadingTokenPtr);
			CleanAfterError();
			}

		/*  Read scale  */
		ScaleFactor = dblstrf (&InputString);

		/*  Determine which kind of potential is being scale	*/
		LeadingTokenPtr = strhed (&InputString);
		__parse_pot (&PotType, &AtomType, LeadingTokenPtr, &InputString);

		/*
		Get pointer to correct potential data
		(Invalid results handled inside routine)
		*/
		e = GetPotPtr (PotType, AtomType);

		/*  Scale input  */
		if (ScaleInput)
			{
			e->x0 *= ScaleFactor;
			e->xn *= ScaleFactor;
			e->dx /= ScaleFactor;
			}

		/*  Scale output	*/
		else
			{
			LOOP (itable, 6*e->n)
				{
				e->a[itable] *= ScaleFactor;
				}
			}

		return;
		}


	/*  Clear Potential	*/
	else if (!strcmpi(LeadingTokenPtr,"CLEAR"))
		{
		/*  Determine which kind of potential is being scale	*/
		LeadingTokenPtr = strhed (&InputString);
		__parse_pot (&PotType, &AtomType, LeadingTokenPtr, &InputString);

		/*
		Select corret potential structure
		(Invalid results handled inside routine)
		*/
		e = GetPotPtr (PotType, AtomType);

		/*  Set variables to zero	*/
		e->x0 = 0.0;
		e->xn = 0.0;
		e->dx = 0.0;
		e->n	= 0;

		/*  else free storage  */
		FREE (e->a)

		return;
		}


	/*  Write Potential	*/
	else if (!strcmpi(LeadingTokenPtr,"WRITE"))
		{

		/*  Call routine for writing potential  */
		WritePotential (InputString);

		/*  Return to calling program  */
		return;
		}


	/*  "Default" potential command (read potential data)  */

	/*  TEST FOR "FORMULA" OPTION  */
	FormulaOptionChosen = ( !strcmpi(LeadingTokenPtr, "FORMULA") );
	if (FormulaOptionChosen)
		{

		/*  Get next token  */
		LeadingTokenPtr = strhed (&InputString);

		/*
		Test if next token is number, this indicates number of strings
		in formula.  A value of zero indicates no number which means a default
		of one.
		*/
		NumFormula = atoi (LeadingTokenPtr);

		/*  No number - use default  */
		if (NumFormula==0)
			NumFormula = 1;

		/*  .. else number present, strip it off	*/
		else
			LeadingTokenPtr = strhed (&InputString);

		}



	/*  "STANDARD"  POTENTIAL COMMAND  */
	__parse_pot (&PotType, &AtomType, LeadingTokenPtr, &InputString);
	e = GetPotPtr (PotType, AtomType);

	/*  FREE PREVIOUS EAMDATA_T  */
	if (e!=NULL && e->a!=NULL)
		{
		FREE (e->a)
		}

	/*  PARSE OPTIONS  */
	e->n	= intstrf (&InputString);
	e->x0 = dblstrf (&InputString);
	e->xn = dblstrf (&InputString);

	/*  SCALE DISTANCE IF PAIR OR DENS	*/
	if (PotType==POT_PAIR || PotType==POT_DENS)
		{
		e->x0 *= 1e-8;
		e->xn *= 1e-8;
		}
	e->dx = (e->n - 1) / (e->xn - e->x0);
	n		= e->n;

	/*  Allocate room for temporary and permanant potential table storage  */
	ALLOCATE (atmp, double, n+4)
	ALLOCATE (e->a, double, 6*n)

	/* Make it so array a starts at index -2	*/
	a		= atmp+2;

	/*  Fill table from formula  */
	if (FormulaOptionChosen)
		{
		double InputScale;
		double ScaleFactor;
		double IndepedentVariable;
		char	 VariableStringBuffer[NUM_BUFFER];
		char	 *StringPtr;
		int	 ErrorCode;
		char	 **FormulaBuffer=NULL;

		/*  Allocate formula storage	*/
		ALLOCATE (FormulaBuffer, char *, NumFormula)
		LOOP (iformula, NumFormula)
			ALLOCATE (FormulaBuffer[iformula], char, NUM_BUFFER)

		/*  Read formula string from input	*/
		LOOP (iformula, NumFormula)
			m_gets_list(FormulaBuffer[iformula],NUM_BUFFER,inlist);

		/*  Calculate scale factor for obtaining independent variable	*/
		ScaleFactor = 1.0/e->dx;

		/*
		Corrent input units from cm (units of e->xn and e->x0)
		to Angstroms (units for formula if using Pair or Dens)
		(This way user can write formula in terms of Angstroms
		instead of cm)
		*/
		if (PotType==POT_PAIR || PotType==POT_DENS)
			InputScale = 1e8;
		else
			InputScale = 1.0;

		/*  Loop over grid values of independent variable	*/
		LOOP (itable, e->n)
			{
			/*  Write INDEPEDENT VARIABLE to parsing variable R  */
			IndepedentVariable = InputScale * (e->x0 + ScaleFactor*itable);


			/*  Form string R=nnn for passage to dblstrf()	*/
			sprintf
				(VariableStringBuffer, "R=%le", IndepedentVariable);

			/*
			Call to dblstrf() evaluates expression
			and sets value of variable R internally
			*/
			StringPtr = VariableStringBuffer;
			dblstrf (&StringPtr);

			/*
			Evaluate formula
			(will use value of R passed via dblstrf()
			*/
			iformula = 0;
			ErrorCode	= FALSE;
			while (iformula<NumFormula && !ErrorCode)
				{
				StringPtr = FormulaBuffer[iformula];
				a[itable] = evalform (&StringPtr, &ErrorCode);
				iformula++;
				}

			/*  Error checking  */
			if (ErrorCode)
				{
				printf ("ERROR (em_read_potential):  Cannot parse formula\n");
				printf ("   Parser returns following message: %s\n",
					parsemsg(ErrorCode));
				CleanAfterError();
				}
			}

		/*  Free formula buffer  */
		LOOP (iformula, NumFormula)
			FREE (FormulaBuffer[iformula])
		FREE (FormulaBuffer)

		}