xmd Code

/*
xmd - molecular dynamics for metals and ceramics
By Jonathan Rifkin <jon.rifkin@uconn.edu>
Copyright 1995-2004 Jonathan Rifkin

This program is free software; you can redistribute it and/or
modify it under the terms of the GNU General Public License
as published by the Free Software Foundation; either version 2
of the License, or (at your option) any later version.

This program is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
GNU General Public License for more details.

You should have received a copy of the GNU General Public License
along with this program; if not, write to the Free Software
Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA  02111-1307, USA.
*/


/*
************************************************************************
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;
extern Simulation_t *s;



/*
************************************************************************
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] +
         ScaledPosition[Z]*ScaledPosition[Z];

      /*  Skip if inside cavity  */
      if (Measure <= 1.0)
         continue;


      /*  Calculate distance to the ellipse center  */
      Radius = 
         sqrt 
         (
         Position[X]*Position[X] + 
         Position[Y]*Position[Y] +
         Position[Z]*Position[Z]
         );

      /*  
      Calculate the penetration into the wall
         - NOTE: Penetration defined as 
              The distance from particle to the ellipsoid wall
              along the direction running to the ellipsoid center.
           NOT
              The distance from particle to closest point on
              ellipsoid wall.
      */
      Penetration = Radius * (1.0 - 1.0/sqrt(Measure));
      

      /*  
      Calculate normal to ellipsoid wall at intersection 
      of wall and a line joining particle and ellipsoid center
      */
      MagGradientSqr = 0.0;
      LOOP (idir, NDIR)
         {
         Gradient[idir]  = InvAxis [idir] * ScaledPosition[idir];
         MagGradientSqr += Gradient[idir] * Gradient[idir];
         }
      
		/*
		Penetration determines force magnitude
      Gradient    determines force direction
		*/
      Force = -a->CavitySpring * Penetration / sqrt(MagGradientSqr);
      TotalForce[ipart][X] += Force * Gradient[X];
      TotalForce[ipart][Y] += Force * Gradient[Y];
      TotalForce[ipart][Z] += Force * Gradient[Z];


#ifdef SPHERE
      /*  Skip if inside cavity  */
      if (Radius<= a->CavityRadius)
         continue;

      /*  Calculate penetration into the cavity wall  */
      Penetration = Radius - a->CavityRadius;

      /*  Calculate force normalized by Radius  */
      Force = -a->CavitySpring * Penetration / Radius;
      TotalForce[ipart][X] += Force * Position[X];
      TotalForce[ipart][Y] += Force * Position[Y];
      TotalForce[ipart][Z] += Force * Position[Z];
#endif

      
      }  /*  LOOP(ipart,a->np)  */
   }  


/*  APPLY EXTERNAL FORCE  */
void ApplyExternalForce (Particle_t *a)
   {
   double (*Coord     )[NDIR];
   double (*TotalForce)[NDIR];
   double *ExternalForce;
   double *ForceOrigin;
   double Displacement;
   double PotentialEnergy;
   int    ipart;
   int    idir;

   /*  Return if no external forces  */
   if (a->ForcePtrList == NULL)
      return;

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

   /*  Loop over all particles subject external forces  */
   PotentialEnergy = 0.0;
   LOOP (ipart, a->np)
   if (a->ForcePtrList[ipart] != NULL)
      {

      /*  Store pointer to External force and origin   */
      ExternalForce = a->ForcePtrList[ipart]->Vector;
      ForceOrigin   = a->ForcePtrList[ipart]->Origin;

      /*  Loop over all components   */
      LOOP (idir, NDIR)
         {
         TotalForce[ipart][idir] += ExternalForce[idir];
         Displacement             = Coord[ipart][idir] - ForceOrigin[idir];
         PotentialEnergy         -= ExternalForce[idir] * Displacement;
         }

      } /*   if (a->ForcePtrList[ipart] != NULL)  */

   /*  Sum contributions to external energy   */
   a->epot += PotentialEnergy;
   }

/*  Apply external springs  */
void ApplyExternalSpring (Particle_t *a)
   {
   double (*Coord)[NDIR];
   double (*Force)[NDIR];
   double *SpringConstant;
   double PotentialEnergy;
   double SpringMag;
   double Dot;
   int ipart;
   int idir;

   /*  Return if no external springs  */
   if (a->SpringPtrList == NULL)
      return;

   /*  Get pointer to coordinates  */
   Coord = (double (*)[NDIR]) a->cur;
   Force = (double (*)[NDIR]) a->f;


   /*  Initialize potential energy   */
   PotentialEnergy = 0.0;

   /*  Loop over all particles  */
   LOOP (ipart, a->np)

   /*  If spring constants assigned to current particle   */
   if (a->SpringPtrList[ipart] != NULL)
      {

      /*  Calculate force components  */
      LOOP (idir, NDIR)
         {

         /*  Store spring constant pointer */
         SpringConstant = a->SpringPtrList[ipart]->Vector;

         /*  Calculate magnitude squared of spring constant  */
         SpringMag =
            (
            SQR(SpringConstant[X]) +
            SQR(SpringConstant[Y]) +
            SQR(SpringConstant[Z])
            );

         /*  Test for zero spring  */
         if (SpringMag==0)
            continue;

         /*  Get magnitude  */
         SpringMag    = sqrt(SpringMag);

         /*  Calculate Dot Product   */
         Dot = 0.0;
         LOOP (idir, NDIR)
            {
            Dot +=
               SpringConstant[idir] *
               (Coord[ipart][idir] - a->SpringPtrList[ipart]->Origin[idir]);
            }

         /*  Scale dot product by magnitude of spring vector  */
         Dot /= SpringMag;

         /*  Calculate force and Potential */
         PotentialEnergy += SQR(Dot);
         LOOP (idir, NDIR)
            Force[ipart][idir] -= SpringConstant[idir] * Dot;

         }


      }

      /*  Add current potential energy to total  */
      a->epot += 0.5 * SpringMag * PotentialEnergy;
   }



void ApplyDamping (Particle_t *a)
	{
	int ipart;
	int idir;
	double Coeff;
	double (*Force   )[NDIR] = (double (*)[NDIR]) a->f;
	double (*Velocity)[NDIR] = (double (*)[NDIR]) a->v;


	/*  Check if damping is one  */
	ASSERT(!(a->UseDamp && a->Damp==NULL))
	if (a->UseDamp==FALSE)
		return;

	/*  Apply damping  */
	LOOP (ipart, a->np)
		{
		Coeff = a->Damp[ipart] / s->dtime;
		LOOP (idir,  NDIR)
			{
			Force[ipart][idir] -= Coeff * Velocity[ipart][idir];
			}
		}

	}



/*
************************************************************************
Local functions
************************************************************************
*/
void ReadCavity (char *InputStr)
   {
   char *TokenStr;
   
   TokenStr = strhed (&InputStr);
   
   /*  Read sphere option  */
   if (!strcmpi(TokenStr, "SPHERE"))
      ReadCavityEllipsoid (InputStr);

   /*  Read ellipsoid option  */
   else if (!strcmpi(TokenStr, "ELLIPSOID"))
      ReadCavityEllipsoid (InputStr);

   /*  Unrecognized option  */
   else
      {
      printf ("ERROR:  Unknown option (%s) to CAVITY.\n",
         TokenStr);
      CleanAfterError();
      }
   }


void ReadCavityEllipsoid (char *InputStr)
   {
   char *TokenStr;

   TokenStr = strhed (&InputStr);

   if (!strcmpi(TokenStr,"OFF"))
      {
      a->UseCavity = FALSE;
      return;
      }

   a->UseCavity = TRUE;

   a->CavitySpring    = dblstrf (&TokenStr);

   a->CavityCenter[X] = 1e-8*dblstrf (&InputStr);
   a->CavityCenter[Y] = 1e-8*dblstrf (&InputStr);
   a->CavityCenter[Z] = 1e-8*dblstrf (&InputStr);

   a->CavityAxis[X]   = 1e-8*dblstrf (&InputStr);
   a->CavityAxis[Y]   = 1e-8*dblstrf (&InputStr);
   a->CavityAxis[Z]   = 1e-8*dblstrf (&InputStr);

   /*  Handle Sphere case  */
   if (a->CavityAxis[Y]==0.0)  a->CavityAxis[Y] = a->CavityAxis[X];
   if (a->CavityAxis[Z]==0.0)  a->CavityAxis[Z] = a->CavityAxis[Y];

      
   /* Sanity Check  */
   if 
   (
   (a->CavityAxis[X] <= 0.0) ||
   (a->CavityAxis[Y] <= 0.0) ||
   (a->CavityAxis[Z] <= 0.0) 
   )
      {
      printf ("ERROR:  Cavity axis (%le,%le,%le) cannot be less than zero.\n",
         a->CavityAxis[X], a->CavityAxis[Y], a->CavityAxis[Z]);
      CleanAfterError();
      }
   if (a->CavitySpring <= 0.0)
      {
      printf ("ERROR:  Cavity spring (%le) cannot be less than zero.\n",
         a->CavitySpring);
      CleanAfterError();
      }

   }