cdeam.c

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

	/*  Read Table  */
	else
		{
		nread = 0;
		while (nread<n && NULL!=m_gets_list(InputBuffer,NUM_BUFFER,inlist))
			{
			LeadingTokenPtr = InputBuffer;
			while (*LeadingTokenPtr)
				{
				a[nread] = dblstrf (&LeadingTokenPtr);
				nread++;
				}
			}

		if (nread<n)
			{
			printf ("ERROR (em_read_potential): Potential table too short.\n");
			CleanAfterError();
			}
		}

	/*  CONVERT ENERGY UNITS  */
	if (PotType==POT_PAIR || PotType==POT_EMBED)
		LOOP(itable, e->n)
			a[itable] /= s->eunit;

	/*  FUNCTION DOMAIN IN REAL SPACE (PAIR AND DENS)	*/
	if (PotType==POT_PAIR || PotType==POT_DENS)
		{
		a[ -1] = 2.0*a[  0] - a[1];
		a[ -2] = 2.0*a[ -1] - a[0];
		a[n  ] = 	 a[n-1];
		a[n+1] = 	 a[n-1];
		}

	/*  FUNCTION DOMAIN IN RHO SPACE - EXTEND LINEARLY IN BOTH DIRECTIONS  */
	else
		{
		a[ -1] = 2.0*a[  0] - a[  1];
		a[ -2] = 2.0*a[ -1] - a[  0];
		a[n  ] = 2.0*a[n-1] - a[n-2];
		a[n+1] = 2.0*a[n	] - a[n-1];
		}

	/*   INTERPOLATE DERIVATIVES	*/
	p = e->a;
	LOOP (k, e->n-1)
		{
		i = k-2;
		j = k-1;
		l = k+1;
		m = k+2;
		n = k+3;
		*p++ = ( 							a[k]									  )	;
		*p++ = ( 2*a[i] -16*a[j]			  + 16*a[l] - 2*a[m] 		  )/24;
		*p++ = (  -a[i] +16*a[j] - 30*a[k] + 16*a[l] -	 a[m] 		  )/24;
		*p++ = (-9*a[i] +39*a[j] - 70*a[k] + 66*a[l] -33*a[m] + 7*a[n])/24;
		*p++ = (13*a[i] -64*a[j] +126*a[k] -124*a[l] +61*a[m] -12*a[n])/24;
		*p++ = (-5*a[i] +25*a[j] - 50*a[k] + 50*a[l] -25*a[m] + 5*a[n])/24;
		}

	/*  Store last value on interval  */
   *p = a[e->n-1];

	/*  Free unneeded allocation	*/
	FREE (atmp)

	/*  SET CUTOFF  */
	CalcCutoff();

	}

/*
************************************************************************
Local Subroutines
************************************************************************
*/

void __parse_pot (int *ptype, int *type,
						char *tokstr, char **instr)
	{
	int t1, t2;

	if 		(!strcmpi(tokstr, "PAIR"))
		{
		*ptype = POT_PAIR;
		t1 = intstrf (instr)-1;
		t2 = intstrf (instr)-1;
		*type = COMBINED_TYPE(t1,t2);
		}
	else if (!strcmpi(tokstr, "DENS"))
		{
		*ptype = POT_DENS;
		*type = intstrf (instr)-1;
		}
	else if (!strcmpi(tokstr, "EMBED"))
		{
		*ptype = POT_EMBED;
		*type = intstrf (instr)-1;
		}
	else
		*ptype = POT_NONE;

	/*  Test for corrent potential  */
	if (*ptype==POT_NONE)
		{
		printf ("ERROR(__parse_pot):  Potential type unrecognized.\n");
		CleanAfterError();
		}

	/*  Test for correct type  */
	if (*ptype==POT_PAIR)
		{
		if (t1>=NumType_m)
			{
			ERROR_PREFIX
			printf ("First type (%i) is out of range [1..%i]\n",
				t1+1, NumType_m);
			CleanAfterError();
			}
		if (t2>=NumType_m)
			{
			ERROR_PREFIX
			printf ("Second type (%i) is out of range [1..%i]\n",
				t2+1, NumType_m);
			CleanAfterError();
			}
		}

	if (*ptype==POT_EMBED || *ptype==POT_DENS)
		{
		if (*type>=NumType_m)
			{
			ERROR_PREFIX
			printf ("Type (%i) is out of range [1..%i]\n",
				*type+1, NumType_m);
			CleanAfterError();
			}
		}

	}

void __write_potstate_bin	(char *fname)
	{
	int i;
	int npair = NumType_m*(NumType_m+1)/2;
	FILE *fout = fopen (fname, "wb");

	if (fout==NULL)
		{
		printf ("ERROR (__write_potstate_bin ):  Cannot open output file %s.\n",
			fname);
		CleanAfterError();
		}
	fwrite (&NumType_m, sizeof(NumType_m), 1, fout);
	for (i=0;i<npair;i++)
		__write_eamtype_bin	(fout, &Pair_m[i]);
	for (i=0;i<NumType_m;i++)
		__write_eamtype_bin	(fout, &Dens_m[i]);
	for (i=0;i<NumType_m;i++)
		__write_eamtype_bin	(fout, &Embed_m[i]);
	fclose (fout);
	}

void __read_potstate_bin  (char *fname)
	{
	int i;
	int NumPair;
	FILE *fin;

	/*  Open input file	*/
	fin = fopen (fname, "rb");

	if (fin==NULL)
		{
		printf ("ERROR (__read_potstate_bin ):  Cannot open input file %s.\n",
			fname);
		CleanAfterError();
		}

	/*  FREE EXISTING POTENTIALS	*/
	__free_eam();

	/*  READ NUMBER OF TYPES  */
	fread (&NumType_m, sizeof(NumType_m), 1, fin);

	/*  ALLOCATE POTENTIALS  */
	NumPair = NumType_m*(NumType_m+1)/2;
	ALLOCATE (Dens_m,  EAMdata_t, NumType_m);
	ALLOCATE (Embed_m, EAMdata_t, NumType_m);
	ALLOCATE (Pair_m,  EAMdata_t, NumPair	);

	/*  READ POTENTIALS		  */
	for (i=0;i<NumPair;i++)
		__read_eamtype_bin  (fin, &Pair_m[i]);
	for (i=0;i<NumType_m;i++)
		__read_eamtype_bin  (fin, &Dens_m[i]);
	for (i=0;i<NumType_m;i++)
		__read_eamtype_bin  (fin, &Embed_m[i]);

	/*  Close input file  */
	fclose (fin);

	/*  SET CUTOFF  */
	CalcCutoff();
	}

void __write_eamtype_bin  (FILE *fout, EAMdata_t *e)
	{
	fwrite (&(e->x0), sizeof(e->x0), 1, 	fout);
	fwrite (&(e->xn), sizeof(e->xn), 1, 	fout);
	fwrite (&(e->n ), sizeof(e->n ), 1, 	fout);
	if (e->n)
		fwrite (e->a, sizeof(e->a[0]), 6*e->n, fout);
	}

void __read_eamtype_bin  (FILE *fout, EAMdata_t *e)
	{

	fread (&(e->x0), sizeof(e->x0), 1,	  fout);
	fread (&(e->xn), sizeof(e->xn), 1,	  fout);
	fread (&(e->n ), sizeof(e->n ), 1,	  fout);
	if (e->n)
		{
		e->dx = (e->n - 1) / (e->xn - e->x0);
		ALLOCATE (e->a, double, 6*e->n);
		fread (e->a,  sizeof(e->a[0]), 6*e->n, fout);
		}
	}


void __read_potstate_text (char *fname)
	{
	int i,NumPair;
	FILE *fin;
	char bufstr[255];


	fin = fopen (fname, "rt");
	if (fin==NULL)
		{
		printf ("ERROR (__read_text):  Cannot open input file %s.\n",
			fname);
		CleanAfterError();
		}

	/*  FREE EXISTING POTENTIALS	*/
	__free_eam();

	/*  READ NUMBER OF TYPES  */
	fgets (bufstr, 255, fin);
	NumType_m = atoi (bufstr);

	/*  ALLOCATE POTENTIALS  */
	NumPair = NumType_m*(NumType_m+1)/2;
	ALLOCATE (Dens_m,  EAMdata_t, NumType_m);
	ALLOCATE (Embed_m, EAMdata_t, NumType_m);
	ALLOCATE (Pair_m,  EAMdata_t, NumPair	);

	/*  READ POTENTIALS		  */
	for (i=0;i<NumPair;i++)
		__read_eamtype_text (fin, &Pair_m[i]);
	for (i=0;i<NumType_m;i++)
		__read_eamtype_text (fin, &Dens_m[i]);
	for (i=0;i<NumType_m;i++)
		__read_eamtype_text (fin, &Embed_m[i]);
	fclose (fin);

	/*  SET CUTOFF  */
	CalcCutoff();
	}


void __write_potstate_text (char *fname)
	{
	int i,npair;
	FILE *fout = fopen (fname, "wt");

	if (fout==NULL)
		{
		printf ("ERROR (__write_text):  Cannot open output file %s.\n",
			fname);
		CleanAfterError();
		}

	/*  WRITE NUMBER OF TYPES	*/
	fprintf (fout, "%i\n", NumType_m);
	npair = NumType_m*(NumType_m+1)/2;
	/*  WRITE POTENTIALS 		*/
	for (i=0;i<npair;i++)
		__write_eamtype_text (fout, &Pair_m[i]);
	for (i=0;i<NumType_m;i++)
		__write_eamtype_text (fout, &Dens_m[i]);
	for (i=0;i<NumType_m;i++)
		__write_eamtype_text (fout, &Embed_m[i]);
	fclose (fout);
	}

void __read_eamtype_text (FILE *fin, EAMdata_t *e)
	{

	int n;
	char *t, bufstr[255];

	fgets (bufstr, 255, fin);
	t = bufstr;
	e->x0 = dblstrf (&t);
	e->xn = dblstrf (&t);
	e->n	= intstrf (&t);
	if (e->n)
		{
		e->dx = (e->n - 1) / (e->xn - e->x0);

		ALLOCATE (e->a, double, 6*e->n);

		n=0;
		while (n < 6*e->n && NULL!=fgets (bufstr, 255, fin) )
			{
			t = bufstr;
			while (*t)
				e->a[n++] = dblstrf (&t);
			}
		}
	}

void __write_eamtype_text (FILE *fout, EAMdata_t *e)
	{
	int i, nline;

	fprintf (fout, "%20.12e %20.12e %10i\n", e->x0, e->xn, e->n);
	if (e->n)
		{
		nline = 0;
		for (i=0;i<6*e->n;i++)
			{
			if (nline==6)
				{
				nline=0;
				fprintf (fout, "\n");
				}
			fprintf (fout, " %27.19le", e->a[i]);
			nline++;
			}
		fprintf (fout, "\n");
		}
	}


/*  FREE POTENTIALS	*/
void __free_eam (void)
	{
	int i, npair;

	npair = NumType_m*(NumType_m+1)/2;
	if (Pair_m)
		{
		LOOP (i, npair)
			if (Pair_m[i].a)
				FREE (Pair_m[i].a)
		FREE (Pair_m)
		}
	if (Dens_m)
		{
		LOOP (i, NumType_m)
			if (Dens_m[i].a)
				FREE (Dens_m[i].a)
		FREE (Dens_m)
		}
	if (Embed_m)
		{
		LOOP (i, NumType_m)
			if (Embed_m[i].a)
				FREE (Embed_m[i].a)
		FREE (Embed_m)
		}
	FREE(DensRc_m)
	FREE(DensRcSqr_m)
	FREE(Rc_m)
	FREE(RcSqr_m)
	}

static void CalcCutoff(void)
	{
	int tc;
	int t1;
	int t2;
	double cut;


	/*  FIND MAXIMUM CUTOFF FOR EACH PAIR OF TYPES	*/
	s->cutoff = 0.0;

	/*  Find maximum Density cutoff for each pair of types  */
	LOOP (t1, NumType_m)
	LOOP (t2, t1+1)
		{

		/*  Type assigned to t1,t2 pair  */
		tc  = COMBINED_TYPE(t1,t2);

		/*  Maximum density cutoff for pair-of-types  */
		cut = Dens_m[t1].xn;
		if (cut < Dens_m[t2].xn)
			cut = Dens_m[t2].xn;
		DensRc_m[tc] = cut;

		/*  Maximum density/pair-potential cutoff for pair-of-types  */
		if (cut < Pair_m[tc].xn)
			{
			Rc_m[tc] = Pair_m[tc].xn;
			}
		else
			{
			Rc_m[tc] = cut;
			}
		DensRcSqr_m[tc] = DensRc_m[tc] * DensRc_m[tc];
		RcSqr_m[tc]     =     Rc_m[tc] *     Rc_m[tc];
		if (s->cutoff < Rc_m[tc])
			s->cutoff = Rc_m[tc];
		}


	}


#ifndef USE_FN_MACRO
double __pfn (int tt, int d, double r)
	{
	EAMdata_t *e;

	/*  Get pointer to pair potential  */
	e = &(Pair_m[tt]);

	/*  Return zero if potential was initialized  */
	if (e->a==NULL)
		return 0.0;

	/*  Return interpolated value  */
	return __fn (e,d,r);
	}


double __dfn (int t, int d, double r)
	{
	EAMdata_t *e;

	/*  Get pointer to pair potential  */
	e = &(Dens_m[t]);

	/*  Return zero if potential was initialized  */
	if (e->a==NULL)
		return 0.0;

	/*  Return interpolated value  */
	return __fn (e, d, r);
	}
#endif

double __efn (int t, int d, double r)
	{
	EAMdata_t *e;

	/*  Get pointer to pair potential  */
	e = &Embed_m[t];

	/*  Return zero if no potential allocated  */
	if (e->a==NULL)
		return 0.0;

	/*  Return last electron density on interval  */
	if (r == e->xn)
		return e->a[6*(e->n-1)];

	/*  Test for electron density too high  */
	if (r >  e->xn)
		{
		printf ("ERROR:  Electron density (%le) out of range [%le, %le].\n",
			r, e->x0, e->xn);
		CleanAfterError();
		}

	/*  Return interpolated value  */
	return __fn (e, d, r);
	}


/*  Interpolate potential tables  */
double __fn (EAMdata_t *e, int d, double r)
	{
	double q, *p;
	int	 m;

	/*
	If potential not allocated, return 0
	*/
	if (e->a==NULL)
		return 0.0;

	/*
	If r above range, return 0.0
	(Good for dens and pair, embed catches this separately and gives error)
	*/
	if (r>=e->xn)
		return 0.0;

	/*
	If r below range, give function value at low end.
	*/
	if (r<e->x0)
		{
		if (d==0)
			return	  e->a[0];
		else if (d==1)
			return	  e->a[1];
		else if (d==2)
			return 2 * e->a[2];
		else
			return 0.0;
		}

	/*  Find position in interpolation table	*/
	m	= (q = e->dx * (r - e->x0) );
	q -= m;

	/*  Out of range, return zero  */
	if (m>=e->n-1)
		return (0.0);

	/*  Point to m'th item (6 numbers per item)  */
	p = &(e->a[((m << 1) + m) << 1]);

	/*  No derivative  */
	if (d==0)
#ifdef ASM_CODE
		return  __poly (q, 5, p);
#else
		return (p[0] + (p[1] + (p[2] + (p[3] + (p[4] + p[5]*q)*q)*q)*q)*q);
#endif

	/*  First Derivative  */
	if (d==1)
		return ( p[1] + (2*p[2] + (3*p[3] + (4*p[4] + 5*p[5]*q)*q)*q)*q )*e->dx;

	/*  Second Derivative  */
	if (d==2)
		return (2*p[2] + (6*p[3] + (12*p[4] + 20*p[5]*q)*q)*q ) * e->dx * e->dx;

	/*  Higher derivatives - return 0  */
	return 0.0;
	}


#pragma argsused
void em_calcforc	(Particle_t *a, Simulation_t *s )
	{
	unsigned char *type, t1, t2, tt;
	int iproc;
	int i,j,k;
	double xd,yd,zd,xb,yb,zb,r,rsq,forc,tfx,tfy,tfz,t;
	double xbox2,ybox2,zbox2, eembed;
	double rc;
	int xsurf = a->surf[0];
	int ysurf = a->surf[1];
	int zsurf = a->surf[2];
	int    ipart;
	int EndInterval;
	int StartInterval;
	EAMdata_t *e;
	ThreadParam_t *ThreadParam = NULL;
	PTHREAD_T     *Thread      = NULL;
	int           NumProc = 1;


	/*  RECALCULATE NEIGHBORS	*/
   a->CalcUniqueNeighbors = TRUE;
	search_all (a);


	/*  SET UP POINTER TO COORDINATES  */
	eembed= 0.0;
	type	= a->type;
	ALLOCATE (Fdrv_m, double, a->np);
	ALLOCATE (Rho_m,  double, a->np);


	/*  ZERO FORCES  */
	LOOP (ipart, a->np)
		{
		a->f[ipart*NDIR+X] = 0.0;
		a->f[ipart*NDIR+Y] = 0.0;
		a->f[ipart*NDIR+Z] = 0.0;
		}

	/*  Zero Internal Stresses  */
	LOOP (i, NVOIGHT)
		a->Stress[i] = 0.0;

	/* CHECK TYPES 								 */
	LOOP (i, a->np)
		if (type[i]>=NumType_m)
			{
			printf ("ERROR (em_calcforc):  Type for particle %i\n", type[i]+1);
			printf ("                      is out of range [1..%i]\n", NumType_m);
			CleanAfterError();
			}


	/*  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;


	/*  Set module-wide pointer to particle structure  */
	a_m = a;

	/*  Get number of processors  */
	NumProc = GetNumProc();

	/*  Allocate parameter structures  */
	ALLOCATE (ThreadParam, ThreadParam_t, NumProc)

	/*  Get atom limits for each thread  */
	LOOP (iproc, NumProc)
		{

		GetAtomSubset
			(
			NumProc, 
			iproc, 
			&ThreadParam[iproc].First,
			&ThreadParam[iproc].Num
			);
		}


	/*
	Rho calculation
     - call threads or functions depending on if threaded or
       non-threaded version
	*/

#ifdef USE_THREAD_CODE

	/*  Allocate thread variables  */
	ALLOCATE (Thread,      PTHREAD_T,     NumProc)

	pthread_mutex_init (&RhoLock_m, NULL);

	LOOP (iproc, NumProc)
		{

		/*  
		Use thread library to call ThreadCalcRho() 
		from separate threads
		*/
		pthread_create
			( 
			&Thread[iproc],
			NULL,
			(void*) &ThreadCalcRho,
			(void*) &ThreadParam[iproc]
			);
		}

	/*  Wait for threads to finish  */
	LOOP (iproc, NumProc)
		{
		pthread_join(Thread[iproc], NULL);
		}
#endif



#ifdef NOT_USE_THREAD_CODE
	/*  Call ThreadCalcRho() directly  */
	LOOP (iproc, NumProc)
		{
		ThreadCalcRho (&ThreadParam[iproc]);
		}
#endif



	/*  CALCULATE EMBEDDING FUNCTION VALUE AT EACH ATOM  */
	LOOP (i, a->np)
		{

		/*  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_m[i] > e->xn)
			{
			printf 
				("ERROR:  Electron density (%le) out of range [%le, %le]",
				Rho_m[i], e->x0, e->xn);
			printf ("for atom %i (type %i).\n", i+1, type[i]+1);
			CleanAfterError();
			}

		/*  Calculate embedding and its derivative  */
		eembed    += __fn (e, 0, Rho_m[i]);
		Fdrv_m[i]  = __fn (e, 1, Rho_m[i]);
		}


	/*
	Force calculation
     - call threads or functions depending on if threaded or
       non-threaded version
	*/


#ifdef USE_THREAD_CODE
	pthread_mutex_init (&ForceLock_m, NULL);


	/*  Intialize Rho Mutex Lock (this place is as good as any?)  */
	Epair_m = 0.0;
	LOOP (iproc, NumProc)
		{


		/*  
		Use thread library to call ThreadCalcForce() 
		from separate threads
		*/
		pthread_create
			( 
			&Thread[iproc],
			NULL,
			(void*) &ThreadCalcForce,
			(void*) &ThreadParam[iproc]
			);
		}

	/*  Wait for threads to finish  */
	LOOP (iproc, NumProc)
		{
		pthread_join(Thread[iproc], NULL);