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.
*/


/*
************************************************************************
History
************************************************************************
*/
/*
11 Aug 1997
   (1)  in _free_eam() use FREE() macro instead of free().  The macro
        calls free() but then sets pointer value to NULL.  Previously
        the pointer was not being set to NULL, and subsequent calls
        due to a POTENTIAL SET EAM 2 command would try to free() the
        non-NULL pointer.
26 Jan 1998
    (1)  Changed all calls to free() into FREE() macros
 5 Jun 1998
    (1)  Store last value of potential in p[e->n-1]
    (2)  When r=rcut have __efn() return value stored above
14 Sep 1998
    (1)  Add warning messages when reading potential (using
         either standard for FORMULA option) without NTABLE,
         X0 or XN.
    (2)  Correct error in GetAtomSubset() which omitted the
         last atom when using odd number of aatoms.  This error 
         would have no effect when the last atom in array
         was also the last atom in the neighbor list, in which
         case it wouldn't be needed because it was handled
         as a neighbor of previous atoms.
         
*/




/*
************************************************************************
Compiler Switches
************************************************************************
*/
/*  Use PC Assembly code for fabs() and poly() routine  */
#define ASM_CODE
#undef  ASM_CODE

/*  Use special DEBUG file to store diagnostic output  */
#define DEBUG
#undef  DEBUG

/*  Use new thread code (May 12, 1998)  */
#define USE_THREAD_CODE

/*  Substitute __efn __dfn with explicit calls to __fn  */
#define USE_FN_MACRO
#undef  USE_FN_MACRO


/*  
If using INTEL, call FPU initialization at 
beginning of each thread function.  May 12, 1998
*/
#ifdef INTEL
#define INIT_INTEL_FPU asm("finit");
#else
#define INIT_INTEL_FPU
#endif



/*
************************************************************************
File Includes
************************************************************************
*/
#include <stdlib.h>
#include <stdio.h>
#include <math.h>
#include "cdthread.h"
#include "strsub.h"
#include "parse.h"
#include "cdhouse.h"
#include "cdsubs.h"
#include "cdsearch.h"
#include "iomngr.h"


/*
************************************************************************
Type Definitions
************************************************************************
*/

/*  STRUCTURE DEFINITIONS   */
typedef struct
{
    double x0,xn,dx;
    int  n;
    double *a;
} EAMdata_t;


/*  Type for passing info to threads and function */
typedef struct
    {
    int  First;
    int  Num;
    }
    ThreadParam_t;



/*
************************************************************************
Defines
************************************************************************
*/

/*  Constants defining potential types  */
#define POT_NONE    0
#define POT_PAIR    1
#define POT_DENS    2
#define POT_EMBED 3

/*  Size of input string buffer */
#define NUM_BUFFER 256

/*  Number of columns of potential data to a line   */
#define MAX_COL 5


/*
************************************************************************
Macros
************************************************************************
*/
#ifdef USE_FN_MACRO
#define __pfn(TYPE,DERIV,RADIUS)  __fn(&(Pair_m[TYPE]),DERIV,RADIUS)
#define __dfn(TYPE,DERIV,RADIUS)  __fn(&(Dens_m[TYPE]),DERIV,RADIUS)
#endif



/*
************************************************************************
Module-Wide Variables
************************************************************************
*/
static size_t fread_return_m;
static char   *fgets_return_m;

/*
************************************************************************
Local Function Prototypes
************************************************************************
*/
void     __parse_pot (int *, int *, char *, char **);
void     __write_eamtype_bin  (FILE *fout, EAMdata_t *e);
void     __read_eamtype_bin  (FILE *fout, EAMdata_t *e);

#ifndef USE_FN_MACRO
double __pfn (int, int, double);
double __dfn (int, int,       double);
#endif

double __efn (int, int,       double);
double __fn  (EAMdata_t *, int, double);
void     __write_potstate_bin  (char *fname);
void     __read_potstate_bin   (char *fname);
void     __write_eamtype_bin   (FILE *fout, EAMdata_t *e);
void     __read_eamtype_bin   (FILE *fout, EAMdata_t *e);
void     __write_potstate_text (char *fname);
void     __read_potstate_text  (char *fname);
void     __write_eamtype_text  (FILE *fout, EAMdata_t *e);
void     __read_eamtype_text   (FILE *fout, EAMdata_t *e);
void     __free_eam (void);
void     __read_dynamo_element (char *);
void     __read_dynamo_alloy   (char *);
void     __funcstore(int, int, int, double, double, BOOLEAN, double *);

EAMdata_t *GetPotPtr        (int, int);
void          WritePotential    (char *);
int       GetFirstType  (int);
int       GetSecondType (int);

static void  CalcCutoff (void);

static void GetAtomSubset  (int, int, int *, int *);
static void ThreadCalcRho  (void *);
static void ThreadCalcForce(void *);


#ifdef ASM_CODE
double __poly (double, int, double *);
#endif
char *convert_fortran_number_str(char *);

/*
************************************************************************
External Variables
************************************************************************
*/
/*  EXTERNAL VARIABLES  */
extern Simulation_t *s;
extern LIST         *inlist;


/*
************************************************************************
Module Variables
************************************************************************
*/
/*  GLOBAL VARIABLES  */
int  NumType_m=0;
static EAMdata_t *Pair_m    = NULL;
static EAMdata_t *Embed_m   = NULL;
static EAMdata_t *Dens_m    = NULL;
static double   *Rc_m       = NULL;
static double   *RcSqr_m     = NULL;
static double   *DensRc_m    = NULL;
static double   *DensRcSqr_m = NULL;
static double  *Rho_m       = NULL;
static double  *Fdrv_m      = NULL;
static Particle_t  *a_m     = NULL;
static double  Epair_m      = 0.0;

/*  Thread Variables  */
#ifdef USE_THREAD_CODE
static PTHREAD_MUTEX_T RhoLock_m;
static PTHREAD_MUTEX_T ForceLock_m;
static PTHREAD_MUTEX_T StressLock_m;
#endif


/*
************************************************************************
Macros
************************************************************************
*/

/*  Redefine fabs() function    */
#ifdef ASM_CODE
double fabs_new(double);
#define FABS fabs_new
#else
#define FABS fabs
#endif


/*
************************************************************************
Exported Subroutines
************************************************************************
*/

#pragma argsused
void em_energy_list (Particle_t *a, Simulation_t *s, int sflag)
    {
    unsigned char *type, t1, t2, tt;
    int i,j,k;
    double xd,yd,zd,xb,yb,zb,r,rsq,t;
    double xbox2,ybox2,zbox2, etemp, *eatom;
    double *c, *ck, *cj;
    double rc;
    double *rho = NULL;
    int xsurf = a->surf[X];
    int ysurf = a->surf[Y];
    int zsurf = a->surf[Z];
    int EndInterval;
    int StartInterval;
    EAMdata_t *e;


#ifdef DEBUG
    FILE *DebugFile;
    DebugFile = fopen ("Debug.fil", "wt");
    if (DebugFile==NULL)
        {
        printf ("ERROR:  Cannot open debug.fil\n");
        CleanAfterError();
        }
#endif

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

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

    /*  Allocate rho    */
    ALLOCATE (rho, double, a->np);


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

    /* INITIALIZE INDIVIDUAL ATOM ENERGIES  */
    /* CHECK TYPES                               */
    for (i=0;i<a->np;i++)
        {
        a->eatom[i] = 0.0;
        if (type[i]>=NumType_m)
            {
            printf ("ERROR (em_energy_list):\n");
            printf ("  Type for particle %i\n", type[i]+1);
            printf ("  is out of range [1..%i]\n", NumType_m);
            CleanAfterError();
            }
        }


    /* INDIVIDUAL PAIR POTENTIALS  */
    c = a->cur;
    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];
            cj = c + (j << 1) + j;
            ck = c + (k << 1) + k;
            t1 = type[j];
            t2 = type[k];
            tt = COMBINED_TYPE(t1,t2);
            rc = Rc_m[tt];
            /*  dislpacement relative to j  */
            xd = FABS(ck[0] - cj[0]);
            if (!xsurf)
                if (xd>xbox2)
                    xd = xb - xd;
            if (xd<rc)
                {
                yd  = FABS(ck[1] - cj[1]);
                if (!ysurf)
                    if (yd>ybox2)
                        yd = yb - yd;
                if (yd<rc)
                    {
                    zd = FABS(ck[2] - cj[2]);
                    if (!zsurf)
                        if (zd>zbox2)
                            zd = zb - zd;
                    if (zd<rc)
                        {
                        rsq = xd*xd + yd*yd + zd*zd;
                        /*  CALCULATE ENERGY CONTRIBUTIONS FOR CURRENT PAIR  */
                        if (rsq<RcSqr_m[tt])
                            {
                            /*  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);
                                }
                            etemp    = __pfn (tt, 0, r);
#ifdef DEBUG
fprintf (DebugFile, "%i %i %14.4e %14.4e\n", t1, t2, r, etemp);
#endif
                            eatom[j] += etemp;
                            eatom[k] += etemp;
                            }
                        }
                    }
                }
            }
        }

    /*  SUM POTENTIAL AND EMBEDDING ENERGIES    */
    a->epot = 0.0;
    for (i=0;i<a->np;i++) {

        /*  Scale eatom by a half  */
        eatom[i] *= 0.5;

        /*  Calculate embedding function  */
        e = &Embed_m[type[i]];

        /*  No embedding function for this type  */
        if (e->a==NULL) continue;

        /*  Electron Density too high  */
        if (rho[i] > e->xn) {
            printf ("ERROR:  Electron density (%le) out of range [%le, %le]",
                rho[i], e->x0, e->xn);
            printf ("for atom %i (type %i).\n", i+1, type[i]+1);
            CleanAfterError();
         }

        /*  Calculate embedding and its derivative  */
        eatom[i] += __fn(e, 0, rho[i]);
        a->epot += eatom[i];
     }

    /*  FREE ARRAYS  */
    FREE (rho)

#ifdef DEBUG
    fclose (DebugFile);
#endif
    }


#pragma argsused
void em_read_potential (char *LeadingTokenPtr, char *InputString) {
    int PotType;
    int AtomType;
    int ndata;
    double x0,xn;
    char InputBuffer[NUM_BUFFER];
    EAMdata_t *e = NULL;
    double *data  = NULL;
    BOOLEAN ScaleInput;
    double  ScaleFactor;
    int   idata;
    BOOLEAN FormulaOptionChosen;
    int   NumFormula;
    int   iformula;

    /*  POTENTIAL INIT - called internally from XMD, not by user. */
    if (!strcmpi(LeadingTokenPtr,"INIT")) {
        int NumPair;

        /*  FREE POTENTIALS */
        __free_eam();

        /*  READ NEW POTENTIAL */
        NumType_m = intstrf (&InputString);
        if (NumType_m<=0)
            NumType_m = 0;

        /*  TEST FOR NumType_m==0   */
        if (NumType_m==0)
            {
            printf ("ERROR (em_read_potential):  Must specify number of types.\n");
            CleanAfterError();
            }

        /*  ALLOCATE POTENTIALS  */
        NumPair   = NumType_m*(NumType_m+1)/2;
        ALLOCATE (Rc_m,        double,    NumPair);
        ALLOCATE (RcSqr_m,     double,    NumPair);
        ALLOCATE (DensRc_m,    double,    NumPair);
        ALLOCATE (DensRcSqr_m, double,    NumPair);
        ALLOCATE (Dens_m,      EAMdata_t, NumType_m);
        ALLOCATE (Embed_m,     EAMdata_t, NumType_m);
        ALLOCATE (Pair_m,      EAMdata_t, NumPair);
        return;

    /*  POTENTIAL WRITE_STATE_TEXT:  experimental  */
    } else if (!strcmpi(LeadingTokenPtr,"WRITE_STATE_TEXT")) {
        __write_potstate_text (InputString);
        return;


    /*  POTENTIAL READ_STATE_TEXT:  experimental  */
    } else if (!strcmpi(LeadingTokenPtr,"READ_STATE_TEXT")) {
        __read_potstate_text  (InputString);
        return;

    /*  POTENTIAL WRITE_STATE_BIN: experimental  */
    } else if (!strcmpi(LeadingTokenPtr,"WRITE_STATE_BIN")) {
        __write_potstate_bin  (InputString);
        return;

    /*  POTENTIAL READ_STATE_BIN: experimental  */
    } else if (!strcmpi(LeadingTokenPtr,"READ_STATE_BIN")) {
        __read_potstate_bin  (InputString);
        return;

    /*  POTENTIAL READ_DYNAMO_ELEMENT: experimental  */
    /*  Read dynamo element file: funcfl() type  */
    } else if (!strcmpi(LeadingTokenPtr,"READ_DYNAMO_ELEMENT")) {
        __read_dynamo_element (InputString);
        return;

    /*  POTENTIAL READ_DYNAMO_ALLOY: experimental  */
    /*  Read dynamo element file: setfl() type  */
    } else if (!strcmpi(LeadingTokenPtr,"READ_DYNAMO_ALLOY")) {
        __read_dynamo_alloy (InputString);
        return;

    /*   POTENTIAL EVAL { PAIR i j|EMBED i|DENS i} nderiv x  */
    } else if (!strcmpi(LeadingTokenPtr,"EVAL")) {
        int nd;
        double x,f;

        LeadingTokenPtr = strhed (&InputString);
        __parse_pot (&PotType, &AtomType, LeadingTokenPtr, &InputString);
        nd = intstrf (&InputString);
        x   = dblstrf (&InputString);

        /* Evaluate potential and adjust units  */
        if (PotType==POT_PAIR) {
            f = __pfn (AtomType, nd, 1e-8*x) * s->eunit;
		}
        else if (PotType==POT_DENS)
            f = __dfn (AtomType, nd, 1e-8*x);
        else
            f = __efn (AtomType, nd,        x) * s->eunit;
        printf ("        %le\n", f);
        return;

    /*  Scale potential */
    /*  POTENTIAL SCALE {INPUT|OUTPUT} scale {PAIR i j|EMBED i|DENS i}  */
    } else if (!strcmpi(LeadingTokenPtr,"SCALE")) {
        /*  Get option - either "INPUT" or "OUTPUT"  */
        LeadingTokenPtr = strhed (&InputString);

        /*  Determine if scaling function input or output   */
        if    (!strcmpi(LeadingTokenPtr, "INPUT") )
            ScaleInput = TRUE;
        else if (!strcmpi(LeadingTokenPtr, "OUTPUT"))
            ScaleInput = FALSE;
        else {
            printf
                ("ERROR(em_read_potential):  Unrecognized option SCALE %s.\n",
                 LeadingTokenPtr);
            CleanAfterError();
        }

        /*  Read scale  */
        ScaleFactor = dblstrf (&InputString);

        /*  Determine which kind of potential is being scale    */
        LeadingTokenPtr = strhed (&InputString);
        __parse_pot (&PotType, &AtomType, LeadingTokenPtr, &InputString);

        /*
        Get pointer to correct potential data
        (Invalid results handled inside routine)
        */
        e = GetPotPtr (PotType, AtomType);

        /*  Scale input  */
        if (ScaleInput)
            {
            e->x0 *= ScaleFactor;
            e->xn *= ScaleFactor;
            e->dx /= ScaleFactor;
            }

        /*  Scale output    */
        else
            {
            LOOP (idata, 6*e->n)
                {
                e->a[idata] *= ScaleFactor;
                }
            }

        return;


    /*  COMMAND:  POTENTIAL MULTIPLY { PAIR i j | EMBED i | DENS i } */
    /*  Multiply potential by independent variable (usually the radius), inverse of
     *  POTENTIAL DIVIDE.  */
    } else if (!strcmpi(LeadingTokenPtr,"MULTIPLY")) {
        double IndependentVariable;
        int    idata;
		double *table = NULL;

        /*  Determine which kind of potential is being scaled   */
        LeadingTokenPtr = strhed (&InputString);
        __parse_pot (&PotType, &AtomType, LeadingTokenPtr, &InputString);

        /*
        Get pointer to correct potential data
        (Invalid results handled inside routine)
        */
        e = GetPotPtr (PotType, AtomType);

		/*  Make new array to store modified potential table  */
		ALLOCATE (table, double, e->n);
		
        /*  Divide every term: f0(x), f1(x), etc. by radius  */
        LOOP (idata, e->n) 
			table[idata] = e->a[6*idata];

        /*  Multiply function values in table by radius  */
        LOOP (idata, e->n) {
			IndependentVariable = e->x0 + idata/e->dx;   /*  e->dx is really 1/dx  */
			/*  Don't divide by zero  */
			if (IndependentVariable==0) continue;
            /*  Convert IndependentVariable to Angstroms  */
            IndependentVariable /= ANGS;
			table[idata] *= IndependentVariable;
        }

	    __funcstore(PotType, AtomType, e->n, e->x0, e->xn, FALSE, table);

		/*  Free table  */
		FREE(table)

        return;


    /*  COMMAND:  POTENTIAL DIVIDE { PAIR i j | EMBED i | DENS i } */
    /*  Divide potential by independent variable (usually the radius), 
     *  made to normalize DYNAMO potential tables  */
    } else if (!strcmpi(LeadingTokenPtr,"DIVIDE")) {
        double IndependentVariable;
		double *table = NULL;
        int    idata;

        /*  Determine which kind of potential is being scaled   */
        LeadingTokenPtr = strhed (&InputString);
        __parse_pot (&PotType, &AtomType, LeadingTokenPtr, &InputString);

        /*
        Get pointer to correct potential data
        (Invalid results handled inside routine)
        */
        e = GetPotPtr (PotType, AtomType);

		/*  Make new array to store modified potential table  */
		ALLOCATE (table, double, e->n);
		
        /*  Divide every term: f0(x), f1(x), etc. by radius  */
        LOOP (idata, e->n) 
			table[idata] = e->a[6*idata];

        /*  Divide function values in table by radius  */
        LOOP (idata, e->n) {
			IndependentVariable = e->x0 + idata/e->dx;   /*  e->dx is really 1/dx  */
			/*  Don't divide by zero  */
			if (IndependentVariable==0) continue;
            /*  Convert IndependentVariable to Angstroms  */
            IndependentVariable /= ANGS;
			table[idata] /= IndependentVariable;
        }

	    __funcstore(PotType, AtomType, e->n, e->x0, e->xn, FALSE, table);

		/*  Free table  */
		FREE(table)

        return;


    /*  Clear Potential */
    /*  POTENTIAL CLEAR  {PAIR i j|EMBED i|DENS i} */
    } else if (!strcmpi(LeadingTokenPtr,"CLEAR")) {
        /*  Determine which kind of potential is being scale    */
        LeadingTokenPtr = strhed (&InputString);
        __parse_pot (&PotType, &AtomType, LeadingTokenPtr, &InputString);

        /*
        Select corret potential structure
        (Invalid results handled inside routine)
        */
        e = GetPotPtr (PotType, AtomType);

        /*  Set variables to zero   */
        e->x0 = 0.0;
        e->xn = 0.0;
        e->dx = 0.0;
        e->n    = 0;

        /*  else free storage  */
        FREE (e->a)

        return;

    /*  Write Potential in text format  */
    /*  POTENTIAL WRITE  [PLOT] {PAIR i j|EMBED i|DENS i} [n [x0 [xn]]] */
    } else if (!strcmpi(LeadingTokenPtr,"WRITE")) {

        /*  Call routine for writing potential  */
        WritePotential (InputString);

        /*  Return to calling program  */
        return;
    }


    /*  "Default" potential command (read potential data)  */

    /*  TEST FOR "FORMULA" OPTION  */
    /*  POTENTIAL FORMULA [nstring] formstr1 [formstr2 [ .. formstrn ..]]  {PAIR i j|EMBED i|DENS i} ntable x0 xn */
    FormulaOptionChosen = ( !strcmpi(LeadingTokenPtr, "FORMULA") );
    if (FormulaOptionChosen) {

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

        /*
        Test if next token is number, this indicates number of strings
        in formula.  A value of zero indicates no number which means a default
        of one.
        */
        NumFormula = atoi (LeadingTokenPtr);

        /*  No number - use default  */
        if (NumFormula==0)
            NumFormula = 1;

        /*  .. else number present, strip it off    */
        else
            LeadingTokenPtr = strhed (&InputString);
    }

    /*  "STANDARD"  POTENTIAL COMMAND  */
    __parse_pot (&PotType, &AtomType, LeadingTokenPtr, &InputString);

    /*  PARSE OPTIONS  */
    ndata = intstrf (&InputString);
    x0    = dblstrf (&InputString);
    xn    = dblstrf (&InputString);

    /*  Sanity check  */
    if (ndata<3) {
        printf ("ERROR:  Potential table must have at least three points.\n");
        CleanAfterError();
    }
    if (x0>=xn) {
        printf ("ERROR:  Potential function argument x0 must be less than xn\n");
        CleanAfterError();
    }

    /*  Allocate room for temporary data storage  */
    ALLOCATE (data, double, ndata)

    /*  Fill table from formula  */
    if (FormulaOptionChosen) {
        double InputScale;
        double ScaleFactor;
        double IndependentVariable;
        char     VariableStringBuffer[NUM_BUFFER];
        char     *StringPtr;
        int  ErrorCode;
        char     **FormulaBuffer=NULL;

        /*  Allocate formula storage    */
        ALLOCATE (FormulaBuffer, char *, NumFormula)
        LOOP (iformula, NumFormula)
            ALLOCATE (FormulaBuffer[iformula], char, NUM_BUFFER)

        /*  Read formula string from input  */
        LOOP (iformula, NumFormula)
            m_gets_list(FormulaBuffer[iformula],NUM_BUFFER,inlist);

        /*  Calculate scale factor for obtaining independent variable   */
        ScaleFactor = (ndata-1)/(xn-x0);

        /*
        Corrent input units from cm (units of e->xn and e->x0)
        to Angstroms (units for formula if using Pair or Dens)
        (This way user can write formula in terms of Angstroms
        instead of cm)
        */
        if (PotType==POT_PAIR || PotType==POT_DENS)
            InputScale = 1e8;
        else
            InputScale = 1.0;

        /*  Loop over grid values of independent variable   */
        LOOP (idata, ndata) {
            /*  Write INDEPEDENT VARIABLE to parsing variable R  */
            IndependentVariable = InputScale * (x0 + ScaleFactor*idata);

            /*  Form string R=nnn for passage to dblstrf()  */
            sprintf
                (VariableStringBuffer, "R=%le", IndependentVariable);

            /*
            Call to dblstrf() evaluates expression
            and sets value of variable R internally
            */
            StringPtr = VariableStringBuffer;
            dblstrf (&StringPtr);

            /*
            Evaluate formula
            (will use value of R passed via dblstrf()
            */
            iformula = 0;
            ErrorCode   = FALSE;
            while (iformula<NumFormula && !ErrorCode) {
                StringPtr = FormulaBuffer[iformula];
                data[idata] = evalform (&StringPtr, &ErrorCode);
                iformula++;
            }

            /*  Error checking  */
            if (ErrorCode) {
                printf ("ERROR (em_read_potential):  Cannot parse formula\n");
                printf ("   Parser returns following message: %s\n",
                    parsemsg(ErrorCode));
                CleanAfterError();
            }
        }

        /*  Free formula buffer  */
        LOOP (iformula, NumFormula)
            FREE (FormulaBuffer[iformula])
        FREE (FormulaBuffer)

    /*  Read Table  */
    /*  POTENTIAL {PAIR i j|EMBED i|DENS i} ntable x0 xn  */
    } else {
        idata = 0;
        while (idata<ndata && NULL!=m_gets_list(InputBuffer,NUM_BUFFER,inlist)) {
            LeadingTokenPtr = InputBuffer;
            while (*LeadingTokenPtr) {
                data[idata] = dblstrf (&LeadingTokenPtr);
                idata++;
            }
        }

        if (idata<ndata) {
            printf ("ERROR (em_read_potential): Potential table too short.\n");
            CleanAfterError();
        }
    }

    /*  Process table data and store in module-wide EAM tables  */
    __funcstore(PotType, AtomType, ndata , x0, xn, TRUE, data);

    /*  Free data storage  */
    FREE (data)

    /*  SET CUTOFF  */
    CalcCutoff();

}


/*
************************************************************************
Local Subroutines
************************************************************************
*/

void __parse_pot (int *ptype, int *type, char *tokstr, char **instr) {
    int t1, t2;

    if      (!strcmpi(tokstr, "PAIR")) {
        *ptype = POT_PAIR;
        t1 = intstrf (instr)-1;
        t2 = intstrf (instr)-1;
        *type = COMBINED_TYPE(t1,t2);
    } else if (!strcmpi(tokstr, "DENS")) {
        *ptype = POT_DENS;
        *type = intstrf (instr)-1;
	} else if (!strcmpi(tokstr, "EMBED")) {
        *ptype = POT_EMBED;
        *type = intstrf (instr)-1;
	} else
        *ptype = POT_NONE;

    /*  Test for corrent potential  */
    if (*ptype==POT_NONE) {
        printf ("ERROR(__parse_pot):  Potential type (%s) unrecognized.\n",tokstr);
        CleanAfterError();
	}

    /*  Test for correct type  */
    if (*ptype==POT_PAIR) {
        if (t1>=NumType_m || t1<0) {
            ERROR_PREFIX
            printf ("First type (%i) is out of range [1..%i]\n", t1+1, NumType_m);
            CleanAfterError();
		}
        if (t2>=NumType_m || t2<0) {
            ERROR_PREFIX
            printf ("Second type (%i) is out of range [1..%i]\n",
                t2+1, NumType_m);
            CleanAfterError();
		}
	}

    if (*ptype==POT_EMBED || *ptype==POT_DENS) {
        if (*type>=NumType_m || *type<0) {
            ERROR_PREFIX
            printf ("Type (%i) is out of range [1..%i]\n",
                *type+1, NumType_m);
            CleanAfterError();
		}
	}

}


void __write_potstate_bin   (char *fname)
    {
    int i;
    int npair = NumType_m*(NumType_m+1)/2;
    FILE *fout = fopen (fname, "wb");

    if (fout==NULL)
        {
        printf ("ERROR (__write_potstate_bin ):  Cannot open output file %s.\n",
            fname);
        CleanAfterError();
        }
    fwrite (&NumType_m, sizeof(NumType_m), 1, fout);
    for (i=0;i<npair;i++)
        __write_eamtype_bin (fout, &Pair_m[i]);
    for (i=0;i<NumType_m;i++)
        __write_eamtype_bin (fout, &Dens_m[i]);
    for (i=0;i<NumType_m;i++)
        __write_eamtype_bin (fout, &Embed_m[i]);
    fclose (fout);
    }

void __read_potstate_bin  (char *fname)
    {
    int i;
    int NumPair;
    FILE *fin;

    /*  Open input file */
    fin = fopen (fname, "rb");

    if (fin==NULL)
        {
        printf ("ERROR (__read_potstate_bin ):  Cannot open input file %s.\n",
            fname);
        CleanAfterError();
        }

    /*  FREE EXISTING POTENTIALS    */
    __free_eam();

    /*  READ NUMBER OF TYPES  */
    fread_return_m = fread (&NumType_m, sizeof(NumType_m), 1, fin);

    /*  ALLOCATE POTENTIALS  */
    NumPair = NumType_m*(NumType_m+1)/2;
    ALLOCATE (Dens_m,  EAMdata_t, NumType_m);
    ALLOCATE (Embed_m, EAMdata_t, NumType_m);
    ALLOCATE (Pair_m,  EAMdata_t, NumPair   );

    /*  READ POTENTIALS       */
    for (i=0;i<NumPair;i++)
        __read_eamtype_bin  (fin, &Pair_m[i]);
    for (i=0;i<NumType_m;i++)
        __read_eamtype_bin  (fin, &Dens_m[i]);
    for (i=0;i<NumType_m;i++)
        __read_eamtype_bin  (fin, &Embed_m[i]);

    /*  Close input file  */
    fclose (fin);

    /*  SET CUTOFF  */
    CalcCutoff();
    }

void __write_eamtype_bin  (FILE *fout, EAMdata_t *e)
    {
    fwrite (&(e->x0), sizeof(e->x0), 1,     fout);
    fwrite (&(e->xn), sizeof(e->xn), 1,     fout);
    fwrite (&(e->n ), sizeof(e->n ), 1,     fout);
    if (e->n)
        fwrite (e->a, sizeof(e->a[0]), 6*e->n, fout);
    }

void __read_eamtype_bin  (FILE *fout, EAMdata_t *e)
    {

    fread_return_m = fread (&(e->x0), sizeof(e->x0), 1,    fout);
    fread_return_m = fread (&(e->xn), sizeof(e->xn), 1,    fout);
    fread_return_m = fread (&(e->n ), sizeof(e->n ), 1,    fout);
    if (e->n)
        {
        e->dx = (e->n - 1) / (e->xn - e->x0);
        ALLOCATE (e->a, double, 6*e->n);
        fread_return_m = fread (e->a,  sizeof(e->a[0]), 6*e->n, fout);
        }
    }


void __read_potstate_text (char *fname)
    {
    int i,NumPair;
    FILE *fin;
    char bufstr[255];


    fin = fopen (fname, "rt");
    if (fin==NULL)
        {
        printf ("ERROR (__read_text):  Cannot open input file %s.\n",
            fname);
        CleanAfterError();
        }

    /*  FREE EXISTING POTENTIALS    */
    __free_eam();

    /*  READ NUMBER OF TYPES  */
    fgets_return_m = fgets (bufstr, 255, fin);
    NumType_m = atoi (bufstr);

    /*  ALLOCATE POTENTIALS  */
    NumPair = NumType_m*(NumType_m+1)/2;
    ALLOCATE (Dens_m,  EAMdata_t, NumType_m);
    ALLOCATE (Embed_m, EAMdata_t, NumType_m);
    ALLOCATE (Pair_m,  EAMdata_t, NumPair   );

    /*  READ POTENTIALS       */
    for (i=0;i<NumPair;i++)
        __read_eamtype_text (fin, &Pair_m[i]);
    for (i=0;i<NumType_m;i++)
        __read_eamtype_text (fin, &Dens_m[i]);
    for (i=0;i<NumType_m;i++)
        __read_eamtype_text (fin, &Embed_m[i]);
    fclose (fin);

    /*  SET CUTOFF  */
    CalcCutoff();
    }


void __write_potstate_text (char *fname)
    {
    int i,npair;
    FILE *fout = fopen (fname, "wt");

    if (fout==NULL)
        {
        printf ("ERROR (__write_text):  Cannot open output file %s.\n",
            fname);
        CleanAfterError();
        }

    /*  WRITE NUMBER OF TYPES   */
    fprintf (fout, "%i\n", NumType_m);
    npair = NumType_m*(NumType_m+1)/2;
    /*  WRITE POTENTIALS        */
    for (i=0;i<npair;i++)
        __write_eamtype_text (fout, &Pair_m[i]);
    for (i=0;i<NumType_m;i++)
        __write_eamtype_text (fout, &Dens_m[i]);
    for (i=0;i<NumType_m;i++)
        __write_eamtype_text (fout, &Embed_m[i]);
    fclose (fout);
    }

void __read_eamtype_text (FILE *fin, EAMdata_t *e)
    {

    int n;
    char *t, bufstr[255];

    fgets_return_m = fgets (bufstr, 255, fin);
    t = bufstr;
    e->x0 = dblstrf (&t);
    e->xn = dblstrf (&t);
    e->n    = intstrf (&t);
    if (e->n)
        {
        e->dx = (e->n - 1) / (e->xn - e->x0);

        ALLOCATE (e->a, double, 6*e->n);

        n=0;
        while (n < 6*e->n && NULL!=fgets (bufstr, 255, fin) )
            {
            t = bufstr;
            while (*t)
                e->a[n++] = dblstrf (&t);
            }
        }
    }

void __write_eamtype_text (FILE *fout, EAMdata_t *e)
    {
    int i, nline;

    fprintf (fout, "%20.12e %20.12e %10i\n", e->x0, e->xn, e->n);
    if (e->n)
        {
        nline = 0;
        for (i=0;i<6*e->n;i++)
            {
            if (nline==6)
                {
                nline=0;
                fprintf (fout, "\n");
                }
            fprintf (fout, " %27.19le", e->a[i]);
            nline++;
            }
        fprintf (fout, "\n");
        }
    }


/*  FREE POTENTIALS */
void __free_eam (void)
    {
    int i, npair;

    npair = NumType_m*(NumType_m+1)/2;
    if (Pair_m)
        {
        LOOP (i, npair)
            if (Pair_m[i].a)
                FREE (Pair_m[i].a)
        FREE (Pair_m)
        }
    if (Dens_m)
        {
        LOOP (i, NumType_m)
            if (Dens_m[i].a)
                FREE (Dens_m[i].a)
        FREE (Dens_m)
        }
    if (Embed_m)
        {
        LOOP (i, NumType_m)
            if (Embed_m[i].a)
                FREE (Embed_m[i].a)
        FREE (Embed_m)
        }
    FREE(DensRc_m)
    FREE(DensRcSqr_m)
    FREE(Rc_m)
    FREE(RcSqr_m)
    }

static void CalcCutoff(void)
    {
    int tc;
    int t1;
    int t2;
    double cut;


    /*  FIND MAXIMUM CUTOFF FOR EACH PAIR OF TYPES  */
    s->cutoff = 0.0;

    /*  Find maximum Density cutoff for each pair of types  */
    LOOP (t1, NumType_m)
    LOOP (t2, t1+1)
        {

        /*  Type assigned to t1,t2 pair  */
        tc  = COMBINED_TYPE(t1,t2);

        /*  Maximum density cutoff for pair-of-types  */
        cut = Dens_m[t1].xn;
        if (cut < Dens_m[t2].xn)
            cut = Dens_m[t2].xn;
        DensRc_m[tc] = cut;

        /*  Maximum density/pair-potential cutoff for pair-of-types  */
        if (cut < Pair_m[tc].xn)
            {
            Rc_m[tc] = Pair_m[tc].xn;
            }
        else
            {
            Rc_m[tc] = cut;
            }
        DensRcSqr_m[tc] = DensRc_m[tc] * DensRc_m[tc];
        RcSqr_m[tc]     =     Rc_m[tc] *     Rc_m[tc];
        if (s->cutoff < Rc_m[tc])
            s->cutoff = Rc_m[tc];
        }


    }


#ifndef USE_FN_MACRO
double __pfn (int tt, int d, double r)
    {
    EAMdata_t *e;

    /*  Get pointer to pair potential  */
    e = &(Pair_m[tt]);

    /*  Return zero if potential was initialized  */
    if (e->a==NULL)
        return 0.0;

    /*  Return interpolated value  */
    return __fn (e,d,r);
    }


double __dfn (int t, int d, double r)
    {
    EAMdata_t *e;

    /*  Get pointer to pair potential  */
    e = &(Dens_m[t]);

    /*  Return zero if potential was initialized  */
    if (e->a==NULL)
        return 0.0;

    /*  Return interpolated value  */
    return __fn (e, d, r);
    }
#endif

double __efn (int t, int d, double r)
    {
    EAMdata_t *e;

    /*  Get pointer to pair potential  */
    e = &Embed_m[t];

    /*  Return zero if no potential allocated  */
    if (e->a==NULL)
        return 0.0;

    /*  Return last electron density on interval  */
    if (r == e->xn)
        return e->a[6*(e->n-1)];

    /*  Test for electron density too high  */
    if (r >  e->xn)
        {
        printf ("ERROR:  Electron density (%le) out of range [%le, %le].\n",
            r, e->x0, e->xn);
        CleanAfterError();
        }

    /*  Return interpolated value  */
    return __fn (e, d, r);
    }


/*  Interpolate potential tables  */
double __fn (EAMdata_t *e, int d, double r) {
    double q, *p;
    int  m;

    /* If potential not allocated, return 0 */
    if (e->a==NULL) return 0.0;

    /* If r above range, return 0.0
    (Good for dens and pair, embed catches this separately and gives error) */
    if (r>=e->xn) return 0.0;

    /* If r below range, give function value at low end.  */
    if (r<e->x0) {
        if      (d==0) return     e->a[0];
        else if (d==1) return     e->a[1];
        else if (d==2) return 2 * e->a[2];
        else           return 0.0;
    }

    /*  Find position in interpolation table    */
    m   = (q = e->dx * (r - e->x0) );
    q -= m;

    /*  Out of range, return zero  */
    if (m>=e->n-1) return 0.0;

    /*  Point to m'th item (6 numbers per item)  */
    p = &(e->a[6*m]);

    /*  No derivative  */
    if (d==0)
#ifdef ASM_CODE
        return  __poly (q, 5, p);
#else
        return (p[0] + (p[1] + (p[2] + (p[3] + (p[4] + p[5]*q)*q)*q)*q)*q);
#endif

    /*  First Derivative  */
    if (d==1)
        return ( p[1] + (2*p[2] + (3*p[3] + (4*p[4] + 5*p[5]*q)*q)*q)*q )*e->dx;

    /*  Second Derivative  */
    if (d==2)
        return (2*p[2] + (6*p[3] + (12*p[4] + 20*p[5]*q)*q)*q ) * e->dx * e->dx;

    /*  Higher derivatives - return 0  */
    return 0.0;
}


#pragma argsused
void em_calcforc    (Particle_t *a, Simulation_t *s )
    {
    unsigned char *type, t1, t2, tt;
    int iproc;
    int i,j,k;
    double xd,yd,zd,xb,yb,zb,r,rsq,forc,tfx,tfy,tfz,t;
    double xbox2,ybox2,zbox2, eembed;
    double rc;
    int xsurf = a->surf[0];
    int ysurf = a->surf[1];
    int zsurf = a->surf[2];
    int    ipart;
    int EndInterval;
    int StartInterval;
    EAMdata_t *e;
    ThreadParam_t *ThreadParam = NULL;
    PTHREAD_T     *Thread      = NULL;
    int           NumProc = 1;


    /*  RECALCULATE NEIGHBORS   */
   a->CalcUniqueNeighbors = TRUE;
    search_all (a);


    /*  SET UP POINTER TO COORDINATES  */
    eembed= 0.0;
    type    = a->type;
    ALLOCATE (Fdrv_m, double, a->np);
    ALLOCATE (Rho_m,  double, a->np);


    /*  ZERO FORCES  */
    LOOP (ipart, a->np)
        {
        a->f[ipart*NDIR+X] = 0.0;
        a->f[ipart*NDIR+Y] = 0.0;
        a->f[ipart*NDIR+Z] = 0.0;
        }

    /*  Zero Internal Stresses  */
    LOOP (i, NVOIGHT)
        a->Stress[i] = 0.0;

    /* CHECK TYPES                               */
    LOOP (i, a->np)
        if (type[i]>=NumType_m)
            {
            printf ("ERROR (em_calcforc):  Type for particle %i\n", type[i]+1);
            printf ("                      is out of range [1..%i]\n", NumType_m);
            CleanAfterError();
            }


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


    /*  Set module-wide pointer to particle structure  */
    a_m = a;

    /*  Get number of processors  */
    NumProc = GetNumProc();

    /*  Allocate parameter structures  */
    ALLOCATE (ThreadParam, ThreadParam_t, NumProc)

    /*  Get atom limits for each thread  */
    LOOP (iproc, NumProc)
        {

        GetAtomSubset
            (
            NumProc, 
            iproc, 
            &ThreadParam[iproc].First,
            &ThreadParam[iproc].Num
            );
        }


    /*
    Rho calculation
     - call threads or functions depending on if threaded or
       non-threaded version
    */

#ifdef USE_THREAD_CODE

    /*  Allocate thread variables  */
    ALLOCATE (Thread,      PTHREAD_T,     NumProc)

    pthread_mutex_init (&RhoLock_m, NULL);

    LOOP (iproc, NumProc)
        {

        /*  
        Use thread library to call ThreadCalcRho() 
        from separate threads
        */
        pthread_create
            ( 
            &Thread[iproc],
            NULL,
            (void*) &ThreadCalcRho,
            (void*) &ThreadParam[iproc]
            );
        }

    /*  Wait for threads to finish  */
    LOOP (iproc, NumProc)
        {
        pthread_join(Thread[iproc], NULL);
        }
#endif



#ifndef USE_THREAD_CODE
    /*  Call ThreadCalcRho() directly  */
    LOOP (iproc, NumProc)
        {
        ThreadCalcRho (&ThreadParam[iproc]);
        }
#endif



    /*  CALCULATE EMBEDDING FUNCTION VALUE AT EACH ATOM  */
    LOOP (i, a->np)
        {

        /*  Calculate embedding function  */
        e = &Embed_m[type[i]];

        /*  No embedding function for this type  */
        if (e->a==NULL)
            continue;

        /*  Electron Density too high  */
        if (Rho_m[i] > e->xn)
            {
            printf 
                ("ERROR:  Electron density (%le) out of range [%le, %le]",
                Rho_m[i], e->x0, e->xn);
            printf ("for atom %i (type %i).\n", i+1, type[i]+1);
            CleanAfterError();
            }

        /*  Calculate embedding and its derivative  */
        eembed    += __fn (e, 0, Rho_m[i]);
        Fdrv_m[i]  = __fn (e, 1, Rho_m[i]);
        }


    /*
    Force calculation
     - call threads or functions depending on if threaded or
       non-threaded version
    */


#ifdef USE_THREAD_CODE
    pthread_mutex_init (&ForceLock_m, NULL);


    /*  Intialize Rho Mutex Lock (this place is as good as any?)  */
    Epair_m = 0.0;
    LOOP (iproc, NumProc)
        {


        /*  
        Use thread library to call ThreadCalcForce() 
        from separate threads
        */
        pthread_create
            ( 
            &Thread[iproc],
            NULL,
            (void*) &ThreadCalcForce,
            (void*) &ThreadParam[iproc]
            );
        }

    /*  Wait for threads to finish  */
    LOOP (iproc, NumProc)
        {
        pthread_join(Thread[iproc], NULL);
        }
    
    /*  Free thread variable  */
    FREE (Thread)
#endif



#ifndef 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 ("POTENTIAL PAIR %i %i  %i  %14.4e  %14.4e\n",
            FirstType, SecondType, NumPoints,
            InputFactor*IntervalStart, InputFactor*IntervalEnd);

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

    } else if (PotType==POT_EMBED) {
        printf ("POTENTIAL 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



/*  Read Dynamo formatted eam potential file for single element  */
void __read_dynamo_element (char *instr) {
    int i, ntotal, atomtype;
    char line[NSTR];
    int    nrho, nr; 
	double drho, dr, cutoff, radius;
    char   *ptr     = NULL;
    double *alldata = NULL;
    /*  Read file name of dynamo file  */
    char *fname = strhed (&instr);
    FILE *fin   = fopen (fname, "r");
    if (fin==NULL) {
        printf ("ERROR:  Cannot open dynamo element file (%s)\n", fname);
        CleanAfterError();
    }
    /*  Skip first line which is a comment  */
    fgets_return_m = fgets(line,NSTR,fin);
    /*  Skip second line which holds mass in Dynamo format - but we enter using MASS command.  */
    fgets_return_m = fgets(line,NSTR,fin);
    /*  Get number of points in real space (for rho and pair functions) 
     *  and in electron density space (for embedding function), and cutoff  */
    fgets_return_m = fgets(line,NSTR,fin);
	convert_fortran_number_str(line);
    sscanf(line,"%d %lg %d %lg %lg", &nrho, &drho, &nr, &dr, &cutoff);
    /*  Sanity Check  */
    if (drho==0) {
        printf ("ERROR:  While reading Dynamo Element potential, desity step size is zero.\n");
        CleanAfterError();
    }
    if (dr==0) {
        printf ("ERROR:  While reading Dynamo Element potential, radius step size is zero.\n");
        CleanAfterError();
    }
    /*  Read in all data points  */
    ntotal = nrho+nr+nr;
    ALLOCATE (alldata, double, ntotal);
    /*  Following loop adapted from lammps' function PairEAM::grab()  */
    i = 0;
    while (i < ntotal) {
        fgets_return_m = fgets(line,NSTR,fin);
		convert_fortran_number_str(line);
        ptr = strtok(line," \t\n");
        alldata[i++] = atof(ptr);
        while (ptr = strtok(NULL," \t\n")) 
            alldata[i++] = atof(ptr);
    }
    fclose(fin);
    /*  Divide pair potential (and dens?) by radius in angstroms, skip r=0  */
    for (i=1;i<nr;i++) {
        radius = i*dr;
        /*  Scale pair potential  */
        alldata[nrho+i] /= radius;
        /*  Scale electron density */
        alldata[nrho+nr+i] /= radius;
    }
    /*  The Embedding Fuction and Pair Potentials are read in eV, convert to ergs 
     *  (these two function tables are the first two stored in alldata[])
     */
    for (i=0;i<nrho+nr;i++) alldata[i] /= EV;
    /*  Store embedding function in module-wide variable  */
    atomtype = 0;
    __funcstore(POT_EMBED, atomtype, nrho, 0.0, (nrho-1)*drho, TRUE, alldata        );
    __funcstore(POT_PAIR,  atomtype, nr,   0.0, (nr-1)*dr,     TRUE, alldata+nrho   );
    __funcstore(POT_DENS,  atomtype, nr,   0.0, (nr-1)*dr,     TRUE, alldata+nrho+nr);

    /*  Free unneeded allocation    */
    FREE (alldata)

    /*  SET CUTOFF  */
    CalcCutoff();
}


/*  Read Dynamo formatted eam potential file for an alloy */
void __read_dynamo_alloy   (char *instr) {
}


/*  Store function info in proper global structure  */
void __funcstore(int PotType, int AtomType, int nd, double x0, double xn, int do_unit_scaling, double *data) {
    int itable, i,j,k,l,m,n;
    EAMdata_t *e = NULL;
    double *atmp = NULL;
    double *a    = NULL;
    double *p    = NULL;
WRITEMSG

    e = GetPotPtr (PotType, AtomType);

    /*  Free previous data */
    if (e!=NULL && e->a!=NULL)  {
		FREE (e->a);
	}

    /*  PARSE OPTIONS  */
    e->n  = nd;
    e->x0 = x0;
    e->xn = xn;

    /*  SCALE DISTANCE IF PAIR OR DENS  */
	if (do_unit_scaling) {
		if (PotType==POT_PAIR || PotType==POT_DENS) {
			e->x0 *= 1e-8;
			e->xn *= 1e-8;
		}
	}
    e->dx = (e->n - 1) / (e->xn - e->x0);

    /*  Allocate room for temporary and permanant potential table storage  */
    ALLOCATE (atmp, double, nd+4)
    ALLOCATE (e->a, double, 6*nd)

    /* Make it so array a[-2] points to start of storage  */
    a       = atmp+2;

    /*  Copy data from input to array a[]  */
    for (itable=0;itable<nd;itable++)  a[itable] = data[itable];

    /*  CONVERT ENERGY UNITS  */
	if (do_unit_scaling) 
		if (PotType==POT_PAIR || PotType==POT_EMBED)
			LOOP(itable, e->n)
				a[itable] /= s->eunit;


    /*  FUNCTION DOMAIN IN REAL SPACE (PAIR AND DENS)   */
    if (PotType==POT_PAIR || PotType==POT_DENS) {
        a[ -1]  = 2.0*a[   0] - a[1];
        a[ -2]  = 2.0*a[  -1] - a[0];
        a[nd ]  =     a[nd-1];
        a[nd+1] =     a[nd-1];

    /*  FUNCTION DOMAIN IN RHO SPACE - EXTEND LINEARLY IN BOTH DIRECTIONS  */
    } else {
        a[  -1] = 2.0*a[   0] - a[   1];
        a[  -2] = 2.0*a[  -1] - a[   0];
        a[nd  ] = 2.0*a[nd-1] - a[nd-2];
        a[nd+1] = 2.0*a[nd  ] - a[nd-1];
        }

    /*   INTERPOLATE DERIVATIVES    */
    /*   Interpolation scheme guarantees:
     *   continuity of function:        f_n(u=1) =       f_n+1(u=0)
     *   continuity of derivative: d/dx f_n(u=1) =  d/dx f_n+1(u=0)
     */
    p = e->a;
    LOOP (k, e->n-1) {
        i = k-2;
        j = k-1;
        l = k+1;
        m = k+2;
        n = k+3;
        *p++ = (                      a[k]                            ) ;
        *p++ = ( 2*a[i] -16*a[j]           + 16*a[l] - 2*a[m]         )/24;
        *p++ = (  -a[i] +16*a[j] - 30*a[k] + 16*a[l] -   a[m]         )/24;
        *p++ = (-9*a[i] +39*a[j] - 70*a[k] + 66*a[l] -33*a[m] + 7*a[n])/24;
        *p++ = (13*a[i] -64*a[j] +126*a[k] -124*a[l] +61*a[m] -12*a[n])/24;
        *p++ = (-5*a[i] +25*a[j] - 50*a[k] + 50*a[l] -25*a[m] + 5*a[n])/24;
        }

    /*  Store last value on interval  */
    *p = a[e->n-1];

    /*  Free unneeded allocation    */
    FREE (atmp)
}



/*
************************************************************************
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

#ifndef 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

#ifndef 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

    }


/*  Convert a the D in a fortran exponent to E to be compatible with C  */
char *convert_fortran_number_str(char *str) {
	int i;
	for (i=0;i<strlen(str);i++)
		if (str[i]=='d' || str[i]=='D') 
			str[i] = 'e';
	return str;
}