cdqu.c

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

#include <stdlib.h>
#include <stdio.h>
#include <math.h>
#include "cdhouse.h"
#include "iomngr.h"
#include "strsub.h"
#include "parse.h"
#include "cdsubs.h"
#include "cdsearch.h"

		  /***********************************************/
		  /**  EXTERNAL VARIABLES (DEFINED IN CDSUBS.C	**/
		  /***********************************************/

extern Particle_t   *a,*b;
extern Simulation_t *s;
extern FILE 		  *out;
extern LIST 		  *inlist;


#define MAXTYPE 3
static double __rc [MAXTYPE][MAXTYPE];
static double __rc2[MAXTYPE][MAXTYPE];
static double __rm [MAXTYPE][MAXTYPE];
static double __f  [MAXTYPE][MAXTYPE];
static double __d  [MAXTYPE][MAXTYPE];

#pragma argsused
void qu_energy_list	(Particle_t *a, Simulation_t *s, int sflag)	{
	unsigned char *type, t1, t2;
	int i,j,k;
	double xd,yd,zd,xb,yb,zb,r,cutsq,rsq;
	double xbox2,ybox2,zbox2, etemp, *eatom;
	double *c, *ck, *cj,  u,t;
	int xsurf = a->surf[0];
	int ysurf = a->surf[1];
	int zsurf = a->surf[2];
	int StartInterval;
	int EndInterval;

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

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

	/*  INITIALIZE VALUES  */
	xb 	= a->bcur[0];
	yb 	= a->bcur[1];
	zb 	= a->bcur[2];
	xbox2 = xb - s->cutoff;
	ybox2 = yb - s->cutoff;
	zbox2 = zb - s->cutoff;

	/* INITIALIZE INDIVIDUAL ATOM ENERGIES  */
	for (i=0;i<a->np;i++)
		a->eatom[i] = 0.0;


	/* INDIVIDUAL PAIR POTENTIALS  */
#ifdef OLD_NEIGHBOR
	for (i=0;i<a->ng;i++)
	{
		j	= a->list1[i];
		k	= a->list2[i];
#else
	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];
#endif
			cj = c + 3*j;
			ck = c + 3*k;
			t1 = type[j];
			t2 = type[k];
			xd = fabs(*ck - *cj);			/*  dislpacement relative to j  */
			if (!xsurf)
				if (xd>xbox2)
					xd = xb - xd;
			if ( xd < __rc[t1][t2])
			{
				yd  = fabs(ck[1] - cj[1]);
				if (!ysurf)
					if (fabs(yd)>ybox2)
						yd = yb - yd;
				rsq = xd*xd + yd*yd;
				if (rsq< (cutsq=__rc2[t1][t2]) )
				{
					zd = fabs(ck[2] - cj[2]);
					if (!zsurf)
						if (fabs(zd)>zbox2)
							zd = zb - zd;
					rsq += zd*zd;
					/*  CALCULATE ENERGY CONTRIBUTIONS FOR CURRENT PAIR  */
					if (rsq<cutsq)  {
						/*  CALCULATE ENERGY CONTRIBUTIONS FOR CURRENT PAIR  */
						/*  CONTRIBUTION TO j FROM k	*/
						r	= sqrt(rsq);
						u	= __f[t1][t2]*(r-__rc[t1][t2]);
						t  = u*u-1;
						etemp = __d[t1][t2] * ( t*t-1 );
						eatom[j] += etemp;
						eatom[k] += etemp;
					}
				}
			}
		}
#ifndef OLD_NEIGHBOR
	}
#endif

	a->epot = 0.0;
	for (i=0;i<a->np;i++) {
		eatom[i] *= 0.5;
		a->epot	+= eatom[i];
	}

}


void qu_calcforc	(Particle_t *a, Simulation_t *s)
{
	unsigned char *type, t1, t2;
	int i,j,k;
	double xd,yd,zd,xb,yb,zb,r,rsq, u,u2, cutsq;
	double xbox2,ybox2,zbox2, etemp;
	double *eatom, *f, *fk, *fj, ft, fxt, fyt, fzt;
	double *c,*ck, *cj, t;
	int xsurf = a->surf[0];
	int ysurf = a->surf[1];
	int zsurf = a->surf[2];
	int StartInterval;
	int EndInterval;


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

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

	/*  INITIALIZE VALUES  */
	xb 	= a->bcur[0];
	yb 	= a->bcur[1];
	zb 	= a->bcur[2];
	xbox2 = xb - s->cutoff;
	ybox2 = yb - s->cutoff;
	zbox2 = zb - s->cutoff;

	/* INITIALIZE INDIVIDUAL ATOM ENERGIES  */
	f = a->f;
	for (i=0;i<a->np;i++)  {
		eatom[i] = 0.0;
		*f++	= 0.0;
		*f++	= 0.0;
		*f++	= 0.0;
	}



	/* INDIVIDUAL PAIR POTENTIALS  */
	f = a->f;
	c = a->cur;
#ifdef OLD_NEIGHBOR
	for (i=0;i<a->ng;i++)
	{
		j	= a->list1[i];
		k	= a->list2[i];
#else
	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];
#endif
			cj = c + 3*j;
			ck = c + 3*k;
			t1 = type[j];
			t2 = type[k];
			xd = ck[0] - cj[0];		/*  dislpacement relative to j  */
			if (!xsurf)  {
				if (fabs(xd)>xbox2) {
					if (xd<0)
						xd += xb;
					else
						xd -= xb;
				}
			}
			if (fabs(xd)< __rc[t1][t2]) {
				yd  = ck[1] - cj[1];
				if (!ysurf)  {
					if (fabs(yd)>ybox2) {
						if (yd<0)
							yd += yb;
						else
							yd -= yb;
					}
				}
				rsq = xd*xd + yd*yd;
				if (rsq< (cutsq=__rc2[t1][t2]) )  {
					zd = ck[2] - cj[2];
					if (!zsurf)  {
						if (fabs(zd)>zbox2) {
							if (zd<0)
								zd += zb;
							else
								zd -= zb;
						}
					}
					rsq += zd*zd;
					/*  CALCULATE ENERGY CONTRIBUTIONS FOR CURRENT PAIR  */
					if (rsq<cutsq)  {
						/*  CALCULATE ENERGY CONTRIBUTIONS FOR CURRENT PAIR  */
						/*  CONTRIBUTION TO j FROM k	*/
						r	 = sqrt(rsq);
						u	 = __f[t1][t2]*(r-__rc[t1][t2]);
						u2  = u*u;
						ft  = __d[t1][t2]*__f[t1][t2] * u * (u2 - 1) / r;
						fxt = ft*xd;
						fyt = ft*yd;
						fzt = ft*zd;
						fk  = f + 3*k;
						fk[0] -= fxt;
						fk[1] -= fyt;
						fk[2] -= fzt;
						fj  = f + 3*j;
						fj[0] += fxt;
						fj[1] += fyt;
						fj[2] += fzt;
						t      = u2-1;
						etemp = __d[t1][t2] * (t*t-1); 
						eatom[j] += etemp;
						eatom[k] += etemp;
					}
				}
			}
		}
#ifndef OLD_NEIGHBOR
	}
#endif

	a->epot = 0.0;
	f = a->f;
	for (i=0;i<a->np;i++) {
		eatom[i] *= 0.5;
		*f++		*= 4;
		*f++		*= 4;
		*f++		*= 4;
		a->epot	+= eatom[i];
	}

}


void qu_read_potential (char *tokstr, char *instr)
{
	int		 type1,type2, i,j;
	double	 scale, t;

	if (!strcmpi(tokstr, "REMOVE"))  {
		type1 = intstrf (&instr);
		type2 = intstrf (&instr);
		if (type1<0 || type1>=NTYPE) {
			printf (
				"ERROR:  Type %i out of range [0..%i].\n", type1, NTYPE);
			CleanAfterError();
		}
		if (type2<0 || type2>=NTYPE) {
			printf (
				"ERROR:  Type %i out of range [0..%i].\n", type2, NTYPE);
			CleanAfterError();
		}
		__rc [type1][type2] = __rc [type2][type1] = 0.0;
		__rc2[type1][type2] = __rc2[type2][type1] = 0.0;
		__rm [type1][type2] = __rm [type2][type1] = 0.0;
		__f  [type1][type2] = __f	[type2][type1] = 0.0;
		__d  [type1][type2] = __d	[type2][type1] = 0.0;
	}
	else if (!strcmpi(tokstr, "CLEAR"))  {
		for (i=0;i<NTYPE;i++)
		for (j=0;j<NTYPE;j++) {
			__rm [i][j] = 0.0;
			__rc [i][j] = 0.0;
			__rc2[i][j] = 0.0;
			__d  [i][j] = 0.0;
		}
	}
	else if (!strcmpi(tokstr, "SCALE")) {
		scale = dblstrf (&instr);
		tokstr = strhed (&instr);
		if (!strcmpi(tokstr, "ALL"))  {
			for (i=0;i<NTYPE;i++)
			for (j=0;j<NTYPE;j++)
				__d  [i][j] *= scale;
		}
		else {
			type1 = intstrf (&tokstr);
			type2 = intstrf (&instr);
			if (type1<0 || type1>=NTYPE) {
				printf (
					"ERROR:  Type %i out of range [0..%i].\n", type1, NTYPE);
				CleanAfterError();
			}
			if (type2<0 || type2>=NTYPE) {
				printf (
					"ERROR:  Type %i out of range [0..%i].\n", type2, NTYPE);
				CleanAfterError();
			}
			__d  [type1][type2] *= scale;
			__d  [type2][type1]	= __d  [type1][type2];
		}
	}
	/*  READ QUARTIC PARAMETERS  */
	else	{
		type1 = intstrf (&tokstr)-1;
		type2 = intstrf (&instr)-1;
		if (type1<0 || type1>=NTYPE) {
			printf (
				"ERROR:  Type %i out of range [1..%i].\n", type1+1, NTYPE);
			CleanAfterError();
		}
		if (type2<0 || type2>=NTYPE) {
			printf (
				"ERROR:  Type %i out of range [1..%i].\n", type2+1, NTYPE);
			CleanAfterError();
		}
		__d  [type1][type2] = __d	[type2][type1] = dblstrf (&instr)/s->eunit;
		__rm [type1][type2] = __rm [type2][type1] = ANGS*dblstrf (&instr);
		t														= ANGS*dblstrf (&instr);
		__rc [type1][type2] = __rc [type2][type1] = t;
		__rc2[type1][type2] = __rc2[type2][type1] = t*t;
		if (__rm[type1][type2]==__rc[type1][type2])	{
			printf (
				"ERROR:  rcut (%f) = rmid(%f).\n",
					__rc[type1][type2]*CM, __rm[type1][type2]*CM );
			CleanAfterError();
		}
		__f[type1][type2] = __f[type2][type1] = 1/(	 __rm[type1][type2]
																  - __rc[type1][type2] );
		/*  RESET POTENTIAL CUTOFFS  */
		s->cutoff = __rc[0][0];
		for (i=0;i<NTYPE;i++)
		for (j=0;j<NTYPE;j++)
			if (s->cutoff<__rc[i][j])
				s->cutoff = __rc[i][j];
	}
}