cdeam.c

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

		}
	
	/*  Free thread variable  */
	FREE (Thread)
#endif



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


	/*  TOTAL ENERGY	*/
	a->epot = Epair_m + eembed;

	/*  RELEASE MEMORY  */
	FREE (ThreadParam)
	FREE (Rho_m)
	FREE (Fdrv_m)
	}

/*
Return pointer to EAM potential type depending on Potential and Atom types
*/
EAMdata_t *GetPotPtr (int PotType, int AtomType)
	{
	if 	  (PotType == POT_PAIR) 	 return &(Pair_m [AtomType]);
	else if (PotType == POT_DENS) 	 return &(Dens_m [AtomType]);
	else if (PotType == POT_EMBED)	 return &(Embed_m[AtomType]);

	/*  If have fallen this far, then there is an error  */
	printf ("ERROR(GetPotPtr):  Incorrect PotType (%i) sent.\n",
		PotType);
	CleanAfterError();

	/*
	This statement not necessary for run time, but to avoid compiler error
	*/
	return NULL;
	}


void WritePotential (char *InputString)
	{
	int PotType;
	int AtomType;
	char *LeadingTokenPtr;
	EAMdata_t *PotPtr;
	BOOLEAN WritePlotFormat;
	int	  NumPoints;
	int	  ipoint;
	int	  NumCol;
	double  InputFactor;
	double  OutputFactor;
	double  IntervalStart;
	double  IntervalEnd;
	double  ScaleFactor;
	double  ValueX;
	double  ValueY;
	int	  FirstType;
	int	  SecondType;

	/*  Point to leading token  */
	LeadingTokenPtr = strhed (&InputString);

	/*  Check for plot keyword  */
	if (!strcmpi(LeadingTokenPtr, "PLOT"))
		{

		/*  Set flag for Plot Format (write both x and y)	*/
		WritePlotFormat = TRUE;

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

	else
		WritePlotFormat = FALSE;


	/*  Determine potential type	*/
	__parse_pot (&PotType, &AtomType, LeadingTokenPtr, &InputString);

	/*  Get Potential Pointer	*/
	PotPtr = GetPotPtr (PotType, AtomType);

	/*  Get number of points to write  */
	if (IsBlank(InputString))
		NumPoints = PotPtr->n;
	else
		NumPoints = intstrf (&InputString);

	/*  Get beginning of interval  */
	if (IsBlank(InputString))
		IntervalStart = PotPtr->x0;
	else
		IntervalStart = dblstrf (&InputString);

	/*  Get end of interval  */
	if (IsBlank(InputString))
		IntervalEnd = PotPtr->xn;
	else
		IntervalEnd = dblstrf (&InputString);

	/*  Convert from cm to angs for pair and dens potentials  */
	if (PotType==POT_PAIR || PotType==POT_DENS)
		InputFactor 	= 1e8;
	 else
		InputFactor 	= 1.0;

	/*  Store factor for converting from ergs to output units  */
	if (PotType==POT_PAIR || PotType==POT_EMBED)
		OutputFactor = s->eunit;
	else
		OutputFactor = 1.0;

	/*  Write header	*/
	if (PotType==POT_PAIR)
		{
		FirstType  = GetFirstType	(AtomType) + 1;
		SecondType = GetSecondType (AtomType) + 1;
		printf ("PAIR %i %i  %i  %14.4e  %14.4e\n",
			FirstType, SecondType, NumPoints,
			InputFactor*IntervalStart, InputFactor*IntervalEnd);
		}

	else if (PotType==POT_DENS)
		{
		printf ("DENS %i  %i  %14.4e  %14.4e\n",
			AtomType+1, NumPoints,
			InputFactor*IntervalStart, InputFactor*IntervalEnd);
		}

	else if (PotType==POT_EMBED)
		{
		printf ("EMBED %i  %i  %14.4e  %14.4e\n",
			AtomType+1, NumPoints,
			InputFactor*IntervalStart, InputFactor*IntervalEnd);
		}

	/*  Calculate factor to convert from point index to x  */
	ScaleFactor = (IntervalEnd-IntervalStart)/(NumPoints-1);

	/*
	Print table  (NOTE: table will not include last NumPoints+1 value)
	*/
	NumCol = 0;
	LOOP (ipoint, NumPoints)
		{

		/*  Find x value	*/
		ValueX = IntervalStart + ipoint * ScaleFactor;

		/*  Find function value  */
		if (PotType==POT_PAIR)
			ValueY = __pfn (AtomType, 0, ValueX);
		if (PotType==POT_DENS)
			ValueY = __dfn (AtomType, 0, ValueX);
		if (PotType==POT_EMBED)
			ValueY = __efn (AtomType, 0, ValueX);

		/*  If writing plot format, include x value	*/
		if (WritePlotFormat)
			printf
			(" %13.4e  %13.4e\n", InputFactor*ValueX, OutputFactor*ValueY);

		/*  Else write MAX_COL values to a line  */
		else
			{
			printf (" %13.4e", OutputFactor*ValueY);
			NumCol++;
			if (NumCol==MAX_COL)
				{
				printf ("\n");
				NumCol = 0;
				}
			}
		}

		/*  Print final new line character if needed  */
		if (NumCol>0)
			printf ("\n");
	}


/*  Get first atom type for combined pair atom type  */
int GetFirstType (int CombinedType)
	{
	int FirstType;
	int SecondType;

	SecondType = GetSecondType (CombinedType);
	FirstType  = CombinedType - SecondType*(SecondType+1)/2;

	return FirstType;
	}

/*  
Get first atom type for combined pair atom type  
Formula is derived like this

c(a,b) is combinatin of type a,b where a>b.

c(a,b) = a*(a+1)/2 + b

8*c(a,  0) = (2*a+1)^2
8*c(a+1,0) = (2*a+3)^2

0.5*(sqrt(8*c(a,  0)+1)-1) = a
0.5*(sqrt(8*c(a+1,0)+1)-1) = a+1

So the integer portion of lhs of above equations will give
a, the fractional part will be determined by b.

*/
int GetSecondType (int CombinedType)
	{
	int SecondType;

	SecondType = ( sqrt( 8*CombinedType + 1 ) - 1 ) /2;

	return SecondType;
	}



#ifdef ASM_CODE
#pragma inline
double __poly (double x, int n, double *p)
	{
	asm	 {
		.386
		/*  Load x	*/
		fld qword ptr x;
		/*  Load n	*/
		mov bx, n;
		shl bx, 3;
		/*  Load and multiply coeff  */
		push ds;
		lds si, p;
		fld qword ptr [si+bx];
		}
	loop:
	asm	 {
		.386
		or bx,bx;
		jz  done;
		mul st[0],st[1];
		sub bx, 8;
		fadd qword ptr [si+bx];
		jmp loop;
		}
	done:
	asm	 {
		.386
		pop ds;
		ffree st[1];
		}
	}


double fabs_new (double x)
	{
	asm	 {
		.386
		/*  Load x	*/
		fld qword ptr x;
		fabs;
		}
	}
#endif


/*
************************************************************************
Thread Functions
************************************************************************
*/
static void ThreadCalcRho (void *ptr)
	{
	int i;
	int j;
	int k;
	int t1;
	int t2;
	int tc;
	int StartInterval;
	int EndInterval;
	double rc;
	double xd;
	double yd;
	double zd;
	double xb;
	double yb;
	double zb;
	double xbox2;
	double ybox2;
	double zbox2;
	double r;
	double rsq;
	double t;
	double *rho = NULL;
	ThreadParam_t *ThreadParam = NULL;
	double (*Coord)[NDIR] = (double (*)[NDIR]) a_m->cur;
	int xsurf = a_m->surf[X];
	int ysurf = a_m->surf[Y];
	int zsurf = a_m->surf[Z];

	/*  Initialize Intel FPU stack  */
	INIT_INTEL_FPU

	/*  Get parameters  */
	ThreadParam = (ThreadParam_t *) ptr;

#ifdef USE_THREAD_CODE
	/*  Initialize storage for rho  */
	ALLOCATE (rho, double, a_m->np)
#endif

#ifdef NOT_USE_THREAD_CODE
	rho = Rho_m;
#endif
	
	/*  INITIALIZE VALUES  */
	xb 	= a_m->bcur[X];
	yb 	= a_m->bcur[Y];
	zb 	= a_m->bcur[Z];
	xbox2 = xb - s->cutoff;
	ybox2 = yb - s->cutoff;
	zbox2 = zb - s->cutoff;


	/*  CALCULATE ELECTRON DENSITY AT EACH ATOM	*/
	for 
	(
	j=ThreadParam->First;
	j<ThreadParam->First+ThreadParam->Num; 
	j++
	)
		{
		StartInterval = a_m->NeighborListIndex[j];
		EndInterval   = StartInterval + a_m->NeighborListLength[j];
		for (i=StartInterval; i<EndInterval; i++)
			{
			k	= a_m->NeighborList[i];
			t1 = a_m->type[j];
			t2 = a_m->type[k];
			tc = COMBINED_TYPE(t1,t2);
			rc = DensRc_m[tc];

			/*  Displacement relative to j  */
			xd = FABS(Coord[k][X]-Coord[j][X]);
			if (!xsurf)
				if (xd>xbox2)
					xd = xb - xd;

			if (xd>=rc)
				continue;

			yd = FABS(Coord[k][Y]-Coord[j][Y]);
			if (!ysurf)
				if (yd>ybox2)
					yd = yb - yd;

			if (yd >= rc)
				continue;

			zd = FABS(Coord[k][Z]-Coord[j][Z]);
			if (!zsurf)
				if (zd>zbox2)
					zd = zb - zd;

			if (zd >= rc)
				continue;

			rsq = xd*xd + yd*yd + zd*zd;

			/*  CALCULATE ENERGY CONTRIBUTIONS FOR CURRENT PAIR  */
			if (rsq >= RcSqr_m[tc])
				continue;

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

#ifdef USE_THREAD_CODE
	/*  Sum local value of rho into module-wide array  */
	pthread_mutex_lock (&RhoLock_m);
	LOOP (i, a_m->np)
		{
		Rho_m[i] += rho[i];
		}
	pthread_mutex_unlock (&RhoLock_m);

	FREE(rho)
#endif
	}


/*  
Calculate sequence of atom indices which is gives 
the requested fraction of neighbors from the neighbor list
   NumProc = number of processors
   iproc   = number of current processor (from [0..N-1])
*/
void GetAtomSubset(int NumProc, int iproc, int *FirstAtom, int *NumAtom)
	{
	int LastAtom;

	ASSERT(NumProc>0)
	
	*FirstAtom = ( (float)  iproc   /NumProc) * a_m->np;
	LastAtom   = ( (float) (iproc+1)/NumProc) * a_m->np;
	*NumAtom   = LastAtom - *FirstAtom;
	}




void ThreadCalcForce(void *ptr)
	{
	int i;
	int j;
	int k;
	int tj;
	int tk;
	int tc;
	int StartInterval;
	int EndInterval;
	double rc;
	double xd;
	double yd;
	double zd;
	double xb;
	double yb;
	double zb;
	double xbox2;
	double ybox2;
	double zbox2;
	double r;
	double rsq;
	double forc;
	double tfx;
	double tfy;
	double tfz;
	double q;
	double *p;
	double Pair0;
	double Pair1;
	double Densj1;
	double Densk1;
	int     m;
	ThreadParam_t *ThreadParam = NULL;
	Coord_t *Coord        = (Coord_t *) a_m->cur;
	Coord_t *Force        = NULL;
	double  *Stress       = NULL;
	double  *Epair        = NULL;
	int xsurf = a_m->surf[X];
	int ysurf = a_m->surf[Y];
	int zsurf = a_m->surf[Z];

	/*  Initialize Intel FPU stack  */
	INIT_INTEL_FPU

	/*  Get parameters  */
	ThreadParam = (ThreadParam_t *) ptr;

#ifdef USE_THREAD_CODE
	/*  Allocate force and stress  */
	ALLOCATE (Force, Coord_t, a_m->np)
	ALLOCATE (Stress, double, NVOIGHT)
	ALLOCATE (Epair,  double,       1)
#endif

#ifdef NOT_USE_THREAD_CODE
	/*  Point to force and stress  */
	Force =  (Coord_t *) a_m->f;
	Stress = (double  *) a_m->Stress;
	Epair  = (double  *) &Epair_m;
#endif

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


	/*  CALCULATE FORCES DO TO PAIR AND ELECTRON DENSITY	*/
	for 
	(
	j=ThreadParam->First;
	j<ThreadParam->First+ThreadParam->Num; 
	j++
	)
		{
		StartInterval = a_m->NeighborListIndex[j];
		EndInterval   = StartInterval + a_m->NeighborListLength[j];
		for (i=StartInterval; i<EndInterval; i++)
			{
			k	= a_m->NeighborList[i];
			tj = a_m->type[j];
			tk = a_m->type[k];
			tc = COMBINED_TYPE(tj,tk);
			rc = Rc_m[tc];

			/*  X displacement  */
			xd = Coord[k][X]- Coord[j][X];
			if (!xsurf)
				if 	  (xd> xbox2)
					xd -= xb;
				else if (xd<-xbox2)
					xd += xb;

			if (FABS(xd) >= rc)
				continue;

			/*  Y displacement  */
			yd = Coord[k][Y]- Coord[j][Y];
			if (!ysurf)
				if 	  (yd> ybox2)
					yd -= yb;
				else if (yd<-ybox2)
					yd += yb;

			if (FABS(yd) >= rc)
				continue;

			/*  Z displacement  */
			zd = Coord[k][Z]- Coord[j][Z];
			if (!zsurf)
				if 	  (zd> zbox2)
					zd -= zb;
				else if (zd<-zbox2)
					zd += zb;

			if (FABS(zd) >= rc)
				continue;

			rsq = xd*xd + yd*yd + zd*zd;

			/*  CALCULATE ENERGY CONTRIBUTIONS FOR CURRENT PAIR  */
			if (rsq >= RcSqr_m[tc])
				continue;

			/*  CALCULATE ENERGY CONTRIBUTIONS FOR CURRENT PAIR  */
			r	= sqrt(rsq);

			/*  Calculate Pair Contribution  */

			/*  r beyond cutoff  */
			if (r>=Pair_m[tc].xn)
				{
				Pair0 = 0.0;
				Pair1 = 0.0;
				}
			/*  r before interior cutoff  */
			else if (r<=Pair_m[tc].x0)
				{
				Pair0 = Pair_m[tc].a[0];
				Pair1 = Pair_m[tc].a[1];
				}
			else
				{
				q  = Pair_m[tc].dx * (r - Pair_m[tc].x0);
				m  = (int) q;
				q -= m;
				if (m >= Pair_m[tc].n-1)
					{
					Pair0 = 0.0;
					Pair1 = 0.0;
					}
				else
					{
					p = &(Pair_m[tc].a[6*m]);
					Pair0 =
						(p[0] + (p[1] + (p[2] + (p[3] + (p[4] + p[5]*q)*q)*q)*q)*q);
					Pair1 =
						( p[1] + (2*p[2] + (3*p[3] + (4*p[4] + 5*p[5]*q)*q)*q)*q )*
						Pair_m[tc].dx;
					}
				
				}
				

			/*  Calculate Density Contribution due to tj */
			/*  r beyond cutoff  */
			if (r>=Dens_m[tj].xn)
				{
				Densj1 = 0.0;
				}
			/*  r before interior cutoff  */
			else if (r<=Dens_m[tj].x0)
				{
				Densj1 = Dens_m[tj].a[1];
				}
			else
				{
				q  = Dens_m[tj].dx * (r - Dens_m[tj].x0);
				m  = (int) q;
				q -= m;
				if (m >= Dens_m[tj].n-1)
					{
					Densj1 = 0.0;
					}
				else
					{
					p = &(Dens_m[tj].a[6*m]);
					Densj1 = 
						( p[1] + (2*p[2] + (3*p[3] + (4*p[4] + 5*p[5]*q)*q)*q)*q )*
						Dens_m[tj].dx;
					}
				
				}



			/*  Calculate Density Contribution due to tk */

			/*  If both atoms same type, don't need second calculation */
			if (tj==tk)
				{
				}
			/*  r beyond cutoff  */
			else if (r>=Dens_m[tk].xn)
				{
				Densk1 = 0;
				}
			/*  r before interior cutoff  */
			else if (r<=Dens_m[tk].x0)
				{
				Densk1 = Dens_m[tk].a[1];
				}
			else
				{
				q  = Dens_m[tk].dx * (r - Dens_m[tk].x0);
				m  = (int) q;
				q -= m;
				if (m >= Dens_m[tk].n-1)
					{
					Densk1 = 0.0;
					}
				else
					{
					p = &(Dens_m[tk].a[6*m]);
					Densk1 = 
						( p[1] + (2*p[2] + (3*p[3] + (4*p[4] + 5*p[5]*q)*q)*q)*q )*
						Dens_m[tk].dx;
					}
				
				}

			*Epair  += Pair0;

			if (tj==tk)
				forc = ( Pair1 + (Fdrv_m[j] + Fdrv_m[k])*Densj1 ) / r;
			else
				forc = ( Pair1 + Fdrv_m[j]*Densk1 + Fdrv_m[k]*Densj1 ) / r;

			tfx = forc*xd;
			tfy = forc*yd;
			tfz = forc*zd;
			Force[k][X] -= tfx;
			Force[k][Y] -= tfy;
			Force[k][Z] -= tfz;
			Force[j][X] += tfx;
			Force[j][Y] += tfy;
			Force[j][Z] += tfz;

			/*  Sum internal stresses	 */
			if
			(
			s->StressStepOutput.Save ||
			a_m->BoxMotionAlgorithm!=BMA_NONE
			)
				{

				/*  Only include stress between selected atoms, if applicable  */
				if (
					!s->UseSelectStress || 
					(IS_SELECT2(a_m->Selection,j) && IS_SELECT2(a_m->Selection,k)) 
					)
					{
					Stress[XX] -= tfx*xd;
					Stress[YY] -= tfy*yd;
					Stress[ZZ] -= tfz*zd;
					Stress[YZ] -= tfy*zd;
					Stress[ZX] -= tfz*xd;
					Stress[XY] -= tfx*yd;
					}
				}
			}
		}



#ifdef USE_THREAD_CODE
	pthread_mutex_lock (&ForceLock_m);

	/*  Save results  */
	LOOP (i, a_m->np)
		{
		a_m->f[i*NDIR+X] += Force[i][X];
		a_m->f[i*NDIR+Y] += Force[i][Y];
		a_m->f[i*NDIR+Z] += Force[i][Z];
		}


	a_m->Stress[XX] += Stress[XX];
	a_m->Stress[YY] += Stress[YY];
	a_m->Stress[ZZ] += Stress[ZZ];
	a_m->Stress[YZ] += Stress[YZ];
	a_m->Stress[ZX] += Stress[ZX];
	a_m->Stress[XY] += Stress[XY];

	Epair_m += *Epair;

	pthread_mutex_unlock (&ForceLock_m);

	/*  Free storage  */
	FREE(Epair)
	FREE(Force)
	FREE(Stress)
#endif

	}