/* xmd - molecular dynamics for metals and ceramics By Jonathan Rifkin 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 #include #include #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(); } }