Compute forces for the Embedded Atom Model (EAM). More...

#include "eam.h"
#include <stdlib.h>
#include <string.h>
#include <math.h>
#include <assert.h>
#include "constants.h"
#include "memUtils.h"
#include "parallel.h"
#include "linkCells.h"
#include "CoMDTypes.h"
#include "performanceTimers.h"
#include "haloExchange.h"

Include dependency graph for eam.c:

Go to the source code of this file.

Data Structures
struct	InterpolationObjectSt
	Handles interpolation of tabular data. More...

struct	EamPotentialSt
	Derived struct for an EAM potential. More...

Macros
#define	MAX(A, B) ((A) > (B) ? (A) : (B))

Typedefs
typedef struct InterpolationObjectSt	InterpolationObject
	Handles interpolation of tabular data. More...

typedef struct EamPotentialSt	EamPotential
	Derived struct for an EAM potential. More...

Functions
static int	eamForce (SimFlat *s)
	Calculate potential energy and forces for the EAM potential. More...

static void	eamPrint (FILE file, BasePotential pot)

static void	eamDestroy (BasePotential **pot)

static void	eamBcastPotential (EamPotential *pot)
	Broadcasts an EamPotential from rank 0 to all other ranks. More...

static InterpolationObject *	initInterpolationObject (int n, real_t x0, real_t dx, real_t *data)
	Builds a structure to store interpolation data for a tabular function. More...

static void	destroyInterpolationObject (InterpolationObject **table)

static void	interpolate (InterpolationObject table, real_t r, real_t f, real_t *df)
	Interpolate a table to determine f(r) and its derivative f'(r). More...

static void	bcastInterpolationObject (InterpolationObject **table)
	Broadcasts an InterpolationObject from rank 0 to all other ranks. More...

static void	printTableData (InterpolationObject table, const char fileName)

static void	eamReadSetfl (EamPotential pot, const char dir, const char *potName)
	Reads potential data from a setfl file and populates corresponding members and InterpolationObjects in an EamPotential. More...

static void	eamReadFuncfl (EamPotential pot, const char dir, const char *potName)
	Reads potential data from a funcfl file and populates corresponding members and InterpolationObjects in an EamPotential. More...

static void	fileNotFound (const char callSite, const char filename)

static void	notAlloyReady (const char *callSite)

static void	typeNotSupported (const char callSite, const char type)

BasePotential *	initEamPot (const char dir, const char file, const char *type)
	Allocate and initialize the EAM potential data structure. More...

Detailed Description

Compute forces for the Embedded Atom Model (EAM).

The Embedded Atom Model (EAM) is a widely used model of atomic interactions in simple metals.

http://en.wikipedia.org/wiki/Embedded_atom_model

In the EAM, the total potential energy is written as a sum of a pair potential and the embedding energy, F:

$U = \sum_{ij} \varphi(r_{ij}) + \sum_i F({\bar\rho_i})$

The pair potential $\varphi_{ij}$ is a two-body inter-atomic potential, similar to the Lennard-Jones potential, and $F(\bar\rho)$ is interpreted as the energy required to embed an atom in an electron field with density $\bar\rho$ . The local electon density at site i is calulated by summing the "effective electron density" due to all neighbors of atom i:

$\bar\rho_i = \sum_j \rho_j(r_{ij})$

The force on atom i, ${\bf F}_i$ is given by

$\begin{eqnarray*} {\bf F}_i & = & -\nabla_i \sum_{jk} U(r_{jk})\\ & = & - \sum_j\left\{ \varphi'(r_{ij}) + [F'(\bar\rho_i) + F'(\bar\rho_j)]\rho'(r_{ij}) \right\} \hat{r}_{ij} \end{eqnarray*}$

where primes indicate the derivative of a function with respect to its argument and $\hat{r}_{ij}$ is a unit vector in the direction from atom i to atom j.

The form of this force expression has two significant consequences. First, unlike with a simple pair potential, it is not possible to compute the potential energy and the forces on the atoms in a single loop over the pairs. The terms involving $F'(\bar\rho)$ cannot be calculated until $\bar\rho$ is known, but calculating $\bar\rho$ requires a loop over the pairs. Hence the EAM force routine contains three loops.

Loop over all pairs, compute the two-body interaction and the electron density at each atom
Loop over all atoms, compute the embedding energy and its derivative for each atom
Loop over all pairs, compute the embedding energy contribution to the force and add to the two-body force

The second loop over pairs doubles the data motion requirement relative to a simple pair potential.

The second consequence of the force expression is that computing the forces on all atoms requires additional communication beyond the coordinates of all remote atoms within the cutoff distance. This is again because of the terms involving $F'(\bar\rho_j)$ . If atom j is a remote atom, the local task cannot compute $\bar\rho_j$ . (Such a calculation would require all the neighbors of atom j, some of which can be up to 2 times the cutoff distance away from a local atom—outside the typical halo exchange range.)

To obtain the needed remote density we introduce a second halo exchange after loop number 2 to communicate $F'(\bar\rho)$ for remote atoms. This provides the data we need to complete the third loop, but at the cost of introducing a communication operation in the middle of the force routine.

At least two alternate methods can be used to deal with the remote density problem. One possibility is to extend the halo exchange radius for the atom exchange to twice the potential cutoff distance. This is likely undesirable due to large increase in communication volume. The other possibility is to accumulate partial force terms on the tasks where they can be computed. In this method, tasks will compute force contributions for remote atoms, then communicate the partial forces at the end of the halo exchange. This method has the advantage that the communication is deffered until after the force loops, but the disadvantage that three times as much data needs to be set (three components of the force vector instead of a single scalar $F'(\bar\rho)$ .

Definition in file eam.c.

Macro Definition Documentation

#define MAX	(	A,
		B
	)	((A) > (B) ? (A) : (B))

Definition at line 103 of file eam.c.

Typedef Documentation

typedef struct EamPotentialSt EamPotential

Derived struct for an EAM potential.

Uses table lookups for function evaluation. Polymorphic with BasePotential.

See Also: BasePotential

typedef struct InterpolationObjectSt InterpolationObject

Handles interpolation of tabular data.

See Also: initInterpolationObject; interpolate

Function Documentation

void bcastInterpolationObject ( InterpolationObject ** table )

static

Broadcasts an InterpolationObject from rank 0 to all other ranks.

It is commonly the case that the data needed to create the interpolation table is available on only one task (for example, only one task has read the data from a file). Broadcasting the table eliminates the need to put broadcast code in multiple table readers.

See Also: eamBcastPotential

Definition at line 544 of file eam.c.

Here is the call graph for this function:

Here is the caller graph for this function:

void destroyInterpolationObject ( InterpolationObject ** table )

static

Definition at line 476 of file eam.c.

Here is the call graph for this function:

Here is the caller graph for this function:

void eamBcastPotential ( EamPotential * pot )

static

Broadcasts an EamPotential from rank 0 to all other ranks.

If the table coefficients are read from a file only rank 0 does the read. Hence we need to broadcast the potential to all other ranks.

Definition at line 407 of file eam.c.

Here is the call graph for this function:

Here is the caller graph for this function:

void eamDestroy ( BasePotential ** pot )

static

Definition at line 389 of file eam.c.

Here is the call graph for this function:

Here is the caller graph for this function:

int eamForce ( SimFlat * s )

static

Calculate potential energy and forces for the EAM potential.

Three steps are required:

Loop over all atoms and their neighbors, compute the two-body interaction and the electron density at each atom
Loop over all atoms, compute the embedding energy and its derivative for each atom
Loop over all atoms and their neighbors, compute the embedding energy contribution to the force and add to the two-body force

Definition at line 215 of file eam.c.

Here is the call graph for this function:

Here is the caller graph for this function:

void eamPrint	(	FILE *	file,
		BasePotential *	pot
	)

static

Definition at line 377 of file eam.c.

Here is the caller graph for this function:

void eamReadFuncfl	(	EamPotential *	pot,
		const char *	dir,
		const char *	potName
	)

static

Reads potential data from a funcfl file and populates corresponding members and InterpolationObjects in an EamPotential.

funcfl is a file format for tabulated potential functions used by the original EAM code DYNAMO. A funcfl file contains an EAM potential for a single element.

The contents of a funcfl file are:

Line Num	Description
1	comments
2	elem amass latConstant latType
3	nrho drho nr dr rcutoff
4	embedding function values F(rhobar) starting at rhobar=0
...	(nrho values. Multiple values per line allowed.)
x'	electrostatic interation Z(r) starting at r=0
...	(nr values. Multiple values per line allowed.)
y'	electron density values rho(r) starting at r=0
...	(nr values. Multiple values per line allowed.)

Where:

elem : atomic number for this element
amass : atomic mass for this element in AMU
latConstant : lattice constant for this elemnent in Angstroms
lattticeType : lattice type for this element (e.g. FCC)
nrho : number of values for the embedding function, F(rhobar)
drho : table spacing for rhobar
nr : number of values for Z(r) and rho(r)
dr : table spacing for r in Angstroms
rcutoff : potential cut-off distance in Angstroms

funcfl format stores the "electrostatic interation" Z(r). This needs to be converted to the pair potential phi(r). using the formula

$phi = Z(r) * Z(r) / r$

NB: phi is not defined for r = 0

Z(r) is in atomic units (i.e., sqrt[Hartree * bohr]) so it is necesary to convert to eV.

F(rhobar) is in eV.

Definition at line 752 of file eam.c.

Here is the call graph for this function:

Here is the caller graph for this function:

void eamReadSetfl	(	EamPotential *	pot,
		const char *	dir,
		const char *	potName
	)

static

Reads potential data from a setfl file and populates corresponding members and InterpolationObjects in an EamPotential.

setfl is a file format for tabulated potential functions used by the original EAM code DYNAMO. A setfl file contains EAM potentials for multiple elements.

The contents of a setfl file are:

Line Num	Description
1 - 3	comments
4	ntypes type1 type2 ... typen
5	nrho drho nr dr rcutoff
F, rho	Following line 5 there is a block for each atom type with F, and rho.
b1	ielem(i) amass(i) latConst(i) latType(i)
b2	embedding function values F(rhobar) starting at rhobar=0
...	(nrho values. Multiple values per line allowed.)
bn	electron density, starting at r=0
...	(nr values. Multiple values per line allowed.)
repeat	Return to b1 for each atom type.
phi	phi_ij for (1,1), (2,1), (2,2), (3,1), (3,2), (3,3), (4,1), ...,
p1	pair potential between type i and type j, starting at r=0
...	(nr values. Multiple values per line allowed.)
repeat	Return to p1 for each phi_ij

Where:

ntypes : number of element types in the potential
nrho : number of points the embedding energy F(rhobar)
drho : table spacing for rhobar
nr : number of points for rho(r) and phi(r)
dr : table spacing for r in Angstroms
rcutoff : cut-off distance in Angstroms
ielem(i) : atomic number for element(i)
amass(i) : atomic mass for element(i) in AMU
latConst(i) : lattice constant for element(i) in Angstroms
latType(i) : lattice type for element(i)

setfl format stores r*phi(r), so we need to converted to the pair potential phi(r). In the file, phi(r)*r is in eV*Angstroms. NB: phi is not defined for r = 0

F(rhobar) is in eV.

Definition at line 634 of file eam.c.

Here is the call graph for this function:

Here is the caller graph for this function:

void fileNotFound	(	const char *	callSite,
		const char *	filename
	)

static

Definition at line 819 of file eam.c.

Here is the caller graph for this function:

BasePotential* initEamPot	(	const char *	dir,
		const char *	file,
		const char *	type
	)

read

Allocate and initialize the EAM potential data structure.

Parameters

[in]	dir	The directory in which potential table files are found.
[in]	file	The name of the potential table file.
[in]	type	The file format of the potential file (setfl or funcfl).

Definition at line 171 of file eam.c.

Here is the call graph for this function:

Here is the caller graph for this function:

InterpolationObject * initInterpolationObject	(	int	n,
		real_t	x0,
		real_t	dx,
		real_t *	data
	)

static

Builds a structure to store interpolation data for a tabular function.

Interpolation must be supported on the range $[x_0, x_n]$ , where $x_n = n*dx$ .

See Also: interpolate; bcastInterpolationObject; destroyInterpolationObject

Parameters

[in]	n	number of values in the table.
[in]	x0	minimum ordinate value of the table.
[in]	dx	spacing of the ordinate values.
[in]	data	abscissa values. An array of size n.

Definition at line 452 of file eam.c.

Here is the call graph for this function:

Here is the caller graph for this function:

void interpolate	(	InterpolationObject *	table,
		real_t	r,
		real_t *	f,
		real_t *	df
	)

static

Interpolate a table to determine f(r) and its derivative f'(r).

The forces on the particle are much more sensitive to the derivative of the potential than on the potential itself. It is therefore absolutely essential that the interpolated derivatives are smooth and continuous. This function uses simple quadratic interpolation to find f(r). Since quadric interpolants don't have smooth derivatives, f'(r) is computed using a 4 point finite difference stencil.

Interpolation is used heavily by the EAM force routine so this function is a potential performance hot spot. Feel free to reimplement this function (and initInterpolationObject if necessay) with any higher performing implementation of interpolation, as long as the alternate implmentation that has the required smoothness properties. Cubic splines are one common alternate choice.

Parameters

[in]	table	Interpolation table.
[in]	r	Point where function value is needed.
[out]	f	The interpolated value of f(r).
[out]	df	The interpolated value of df(r)/dr.

Definition at line 512 of file eam.c.

Here is the caller graph for this function:

void notAlloyReady ( const char * callSite )

static

Definition at line 826 of file eam.c.

Here is the caller graph for this function:

void printTableData	(	InterpolationObject *	table,
		const char *	fileName
	)

static

Definition at line 575 of file eam.c.

Here is the call graph for this function:

void typeNotSupported	(	const char *	callSite,
		const char *	type
	)

static

Definition at line 834 of file eam.c.

Here is the caller graph for this function:

Data Structures

Macros

Typedefs

Functions

Detailed Description

Macro Definition Documentation

Typedef Documentation

Function Documentation