cdcons.c

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

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

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



/*
************************************************************************
Local Function Prototypes
************************************************************************
*/
void ReadCavity         (char *);
void ReadCavityEllipsoid(char *);
void ApplyCavityConstraint (Particle_t *a);
void ApplyIndividualConstraint (Particle_t *a);


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


                             /*   READ CONSTRAINT  */

void read_constrain(char *instr)
   {
   int ipart;
	int ctype;
   /*  CONSTRAINT ON FLAG   */
   BOOLEAN cflag = TRUE;
   double xp,yp,zp, xd,yd,zd, mag;
	double Values[6];
	double (*Coord)[NDIR] = (double (*)[NDIR]) a->cur;
   char *tokstr = strhed (&instr);

   /*  CAVITY option  */
   if      (!strcmpi(tokstr,"CAVITY"))
      {
      ReadCavity (instr);
      return;
      }

   /*  OFF   */
   if (!strcmpi(tokstr,"OFF" ))
      cflag = FALSE;

   else
      {

      /* LINE   */
      if      (!strcmpi(tokstr,"LINE"))
         ctype = CONSTR_LINE;
      /*  PLANE  */
      else if (!strcmpi(tokstr,"PLANE"))
         ctype = CONSTR_PLANE;
      else
         {
         printf ("ERROR (read_constrain):  CONSTRAINT TYPE ");
         printf ("(%s) must be LINE, PLANE or OFF.\n", tokstr);
         return;
         }
      }

   /*  IF CONSTRAIN ON   */
   if (cflag)
      {

      /*  Allocate space for array-of-pointer to constraints  */
		if (a->cons==NULL)
	      ALLOCATE (a->cons, Constraint_t*, a->NumPartAlloc)

		LOOP (ipart, a->np)
         {
			if (IS_SELECT(ipart))
            {
				if (a->cons[ipart]==NULL) {
	            ALLOCATE (a->cons[ipart], Constraint_t, 1)
				}

		      /*  Map formula string into Values[] array  */
     			 EvaluateFormula 
					(instr, Coord[ipart], 1e+8, Values, 1, 6);

		      xd = Values[0];
		      yd = Values[1];
		      zd = Values[2];
		      xp = Values[3] * ANGS;
		      yp = Values[4] * ANGS;
		      zp = Values[5] * ANGS;

		      /*  NORMALIZE DIRECTION  */
		      mag = sqrt(xd*xd + yd*yd + zd*zd);
		      if (mag!=0)
		         {
		         xd /= mag;
		         yd /= mag;
		         zd /= mag;
		         }

            a->cons[ipart]->type = ctype;
            a->cons[ipart]->xd   = xd;
            a->cons[ipart]->yd   = yd;
            a->cons[ipart]->zd   = zd;
            a->cons[ipart]->xp   = xp;
            a->cons[ipart]->yp   = yp;
            a->cons[ipart]->zp   = zp;
            }
         }

      /*  APPLY CONSTRAIN NOW  */
      ApplyConstraint (a);
      }

   /*  CONSTRAIN OFF  */
   else if (a->cons)
      {
      int NumCons;

      /*
      Move through list, remove selected constraints,
      count remaining constaints
      */
      NumCons = 0;
		LOOP (ipart, a->np)
         {
         if (a->cons[ipart])
            {

            /*  Free Constaint  */
				if (IS_SELECT(ipart))
					{
               FREE (a->cons[ipart])
					}

            /*  Count constraint  */
            else
               NumCons++;
            }
         }
      /*  IF NO ACTIVE PARTICLE CONSTRAINTS REMOVE a->cons ARRAY   */
      if (NumCons==0)
         FREE (a->cons)
      }

	/*  Re-count  number of degrees of freedom  */
	a->DegFree = GetTotalDOF(a);
   }

/*  Apply family of constraints  */
void ApplyConstraint (Particle_t *a)
   {
   ApplyCavityConstraint (a);
   ApplyIndividualConstraint (a);
   }

/*  APPLY INDIVIDUAL PARTICLE CONSTRAINTS TO EITHER LINE OR PLANE  */
void ApplyIndividualConstraint (Particle_t *a) {
   int i;
   double dot;
	double vdot;
   double fdot;
   Constraint_t *cons;
   double  (*c)[NDIR]  = (double (*)[NDIR]) a->cur;
   double  (*v)[NDIR]  = (double (*)[NDIR]) a->v;
   double  (*f)[NDIR]  = (double (*)[NDIR]) a->f;

   /*  IF CONSTRAINT ARRAY EXISTS  */
   if (a->cons==NULL)
      return;

   LOOP (i, a->np) {
      if (a->cons[i]) {
         cons  = a->cons[i];
         dot   = cons->xd * (c[i][X] - cons->xp);
         dot  += cons->yd * (c[i][Y] - cons->yp);
         dot  += cons->zd * (c[i][Z] - cons->zp);
         if (a->v) 
				vdot = cons->xd*v[i][X] + cons->yd*v[i][Y] + cons->zd*v[i][Z];
			if (a->f) 
				fdot = cons->xd*f[i][X] + cons->yd*f[i][Y] + cons->zd*f[i][Z];
         /*  LINE CONSTRAINT   */
         if (cons->type==CONSTR_LINE) {
            /*  PLACE PARTICLES ON LINES   */
            c[i][X] = dot * cons->xd + cons->xp;
            c[i][Y] = dot * cons->yd + cons->yp;
            c[i][Z] = dot * cons->zd + cons->zp;
				/*  CONSTRAIN VELOCITY PARALLEL TO LINE  */
            if (a->v) {
               v[i][X] = vdot * cons->xd;
               v[i][Y] = vdot * cons->yd;
               v[i][Z] = vdot * cons->zd;
            }
            /*  CONSTRAIN FORCE PARALLEL TO LINE  */
            if (a->f) {
               f[i][X] = fdot * cons->xd;
               f[i][Y] = fdot * cons->yd;
               f[i][Z] = fdot * cons->zd;
            }
         /*  PLANE CONSTRAINT  */
         } else if (cons->type==CONSTR_PLANE) {   /* NOTE: cons[i]->type must be either 0 or 1  */
            /*  DISPLACE PARTICLES PARALLEL TO PLANE NORMAL  */
            c[i][X] -=   dot * cons->xd;
            c[i][Y] -=   dot * cons->yd;
            c[i][Z] -=   dot * cons->zd;
            /*  CONSTRAIN VELOCITY, FORCE PERPENDICULAR TO PLANE NORMAL   */
            if (a->v) {
               v[i][X] -= vdot * cons->xd;
               v[i][Y] -= vdot * cons->yd;
               v[i][Z] -= vdot * cons->zd;
            }
            if (a->f) {
               f[i][X] -= fdot * cons->xd;
               f[i][Y] -= fdot * cons->yd;
               f[i][Z] -= fdot * cons->zd;
            }
         }
      }
   }
}


/*  Apply cavity sphere constraint  */
void ApplyCavityConstraint (Particle_t *a)
   {
   int    ipart;
   int    idir;
   double (*Coord     )[NDIR];
   double (*TotalForce)[NDIR];
   double ScaledPosition[NDIR];
   double Position[NDIR];
   double Gradient[NDIR];
   double InvAxis [NDIR];
	double HalfBox [NDIR];
   double Penetration;
   double Radius;
   double Force;
   double Measure;
   double MagGradientSqr;

   /*  Test for cavity constraint  */
   if (!a->UseCavity)
      return; 

   /*  Set pointer to total force and force origin  */
   Coord      = (double (*)[NDIR]) a->cur;
   TotalForce = (double (*)[NDIR]) a->f;

   /*  Initialize factors  */
   LOOP (idir, NDIR)
      {
      InvAxis[idir] = 1.0 / a->CavityAxis[idir];
		HalfBox[idir] = 0.5 * a->bcur[idir];
      }

   /*  Test all particles  */
   LOOP (ipart, a->np)
      {

      /*  Calculate position relative to cavity center  */
      LOOP (idir, NDIR)
         {
			/*  Get position relative to cavity center  */
         Position[idir] = Coord[ipart][idir] - a->CavityCenter[idir];
			
			/*  Correct for box  */
			if (!a->surf[idir])
				{
				if      (Position[idir] < -HalfBox[idir])
					Position[idir] += a->bcur[idir];
				else if (Position[idir] >  HalfBox[idir])
					Position[idir] -= a->bcur[idir];
				}
	
			/*  Find position in ellipitical coordinates  */
         ScaledPosition[idir] = InvAxis[idir] * Position[idir];
         }

      /*  Calculate measure of position relative to ellipsoid  */
      Measure = 
         ScaledPosition[X]*ScaledPosition[X] + 
         ScaledPosition[Y]*ScaledPosition[Y] +