// EMT.h -- Implements the EMT potential.
// -*- c++ -*-
//
// Copyright (C) 2001-2011 Jakob Schiotz and Center for Individual
// Nanoparticle Functionality, Department of Physics, Technical
// University of Denmark. Email: schiotz@fysik.dtu.dk
//
// This file is part of Asap version 3.
// Asap is released under the GNU Lesser Public License (LGPL) version 3.
// However, the parts of Asap distributed within the OpenKIM project
// (including this file) are also released under the Common Development
// and Distribution License (CDDL) version 1.0.
//
// This program is free software: you can redistribute it and/or
// modify it under the terms of the GNU Lesser General Public License
// version 3 as published by the Free Software Foundation. Permission
// to use other versions of the GNU Lesser General Public License may
// granted by Jakob Schiotz or the head of department of the
// Department of Physics, Technical University of Denmark, as
// described in section 14 of the GNU General Public License.
//
// 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
// and the GNU Lesser Public License along with this program. If not,
// see .
#ifndef _EMT_H
#define _EMT_H
#include "AsapPython.h"
#include "Asap.h"
#include "Atoms.h"
#include "Vec.h"
#include "SymTensor.h"
#include "Potential.h"
#include "EMTParameterProvider.h"
#include "TinyMatrix.h"
#include
using std::vector;
namespace ASAPSPACE {
class NeighborLocator;
struct emt_parameters;
/// The Effective Medium Theory (EMT) potential.
/// The Effective Medium Theory (EMT) as originally implemented in Per
/// Stoltze's ARTwork program, and documented in
/// - K. W. Jacobsen, P. Stoltze and J. K. Norskov,
/// Surf. Sci. vol. 366, p. 394-402 (1996).
class EMT : public Potential
{
public:
/// Create an EMT potential optionally using a parameter provider.
EMT(PyObject *prov);
/// Delete the EMT potential.
virtual ~EMT();
virtual string GetName() const {return "EMT";}
virtual string GetRepresentation() const;
/// Whether we subtract E0 from atomic energies.
///
/// Defines the zero of potential energy; if we don't subtract, zero
/// corresponds to infinite seperation. If we do, it is the bulk fcc
/// crystal at the equilibrium lattice constant.
void SetSubtractE0(bool subtractE0) { this->subtractE0 = subtractE0; }
virtual void SetAtoms(PyObject *pyatoms, Atoms* accessobj = NULL);
virtual const vector &GetForces(PyObject *pyatoms);
virtual const vector &GetPotentialEnergies(PyObject *pyatoms);
virtual double GetPotentialEnergy(PyObject *pyatoms);
virtual const vector &GetStresses(PyObject *pyatoms,
bool virialonly=false);
virtual void GetStress(PyObject *pyatoms, double stress[6],
bool virialonly=false);
/// Check if neighbor lists need an update.
///
/// Return true if the neigbor list needs an update. In parallel
/// simulations, the atoms will cause a communication so it returns
/// true if an update is needed on any processor. In parallel
/// simulations, it also instructs the atoms to communicate their
/// update status.
virtual bool CheckNeighborList();
/// Update the neighbor lists
virtual void UpdateNeighborList();
virtual bool CalcReq_Energy(PyObject *pyatoms);
virtual bool CalcReq_Forces(PyObject *pyatoms);
virtual bool CalcReq_Stress(PyObject *pyatoms);
virtual double GetCutoffRadius() const {return rNbCut;}
virtual double GetLatticeConstant() const;
virtual int GetNumberOfAtoms() const {return nAtoms;}
/// Return a pointer to the EMT "density" sigma1.
const vector > &GetSigma1() {return sigma1;}
/// Return a pointer to the EMT "density" sigma2.
const vector > &GetSigma2() {return sigma2;}
/// Print the EMT parameters.
void PrintParameters();
/// Return the neighbor list.
/// Return a BORROWED reference to the Python object containing
/// neighbor list for the
/// potential, if this type of potential supports it, and if it is
/// defined now. Otherwise, return NULL without setting a Python
/// error.
PyObject *GetNeighborList() const {return nblist_obj;}
/// This potential can be used in parallel simulations
virtual bool Parallelizable() const {return true;}
/// Mark that the atoms have changed the boundary conditions.
virtual void BoundaryConditionsChanged();
/// Print memory usage
virtual long PrintMemory() const;
protected:
/// Initialization of the EMT parameters.
virtual void InitParameters();
/// (Re)allocate storage for forces, energies and intermediate results.
virtual void Allocate();
/// (Re)allocate storage for stresses
virtual void AllocateStress();
/// Calculate type numbers from the atomic numbers.
virtual void CalculateIDs();
/// Calculate sigma1 and perhaps sigma2.
virtual void CalculateSigmas(bool calculatesigma2);
/// Calculate energies
virtual void CalculateEnergies();
/// Calculate energies once sigma1 is known
/// If Epot is NULL, energies are not calculated, only derivatives
/// needed for the force.
virtual void CalculateEnergiesAfterSigmas(bool calcEpot);
/// Calculate the forces
virtual void CalculateForces();
/// Calculate forces in a multicomponent system.
virtual void CalculateForcesAfterEnergies();
/// Calculate forces in a system with only one element.
virtual void CalculateForcesAfterEnergiesSingle();
/// Calculate stresses.
virtual void CalculateStresses(bool virialonly);
/// Allocate the neighbor list.
virtual void CreateNeighborList();
/// sigma_batch does the hard work in CalculateSigmas().
virtual void sigma_batch(int *self, int *other, Vec rnb[],
double *sq_dist, int zs, int zo, int n,
bool calculatesigma2, bool partialupdate = false);
/// force_batch does the hard work in CalculateForcesAfterEnergy().
virtual void force_batch(const int *self, const int *other, const Vec rnb[],
const double sq_dist[], const double dEdss[],
const double dEdso[], int zs, int zo, int n);
void distribute_force(const int *self, const int *other,
const double *df, const Vec *rnb, int n);
protected: // Data
// Atoms *atoms; ///< The atoms we are working on
bool ghostatoms; ///< True if atoms have ghosts.
int nAtoms; ///< The number of (real) atoms.
int nSize; ///< Number of atoms including ghost atoms.
NeighborLocator *nblist; ///< The neighborlist object.
PyObject *nblist_obj;
double driftfactor; ///< Drift factor for the neighbor list.
EMTParameterProvider *provider; ///< The source of the EMT parameters
PyObject *provider_obj; ///< The Python object behind provider.
bool subtractE0; ///< Whether we subtract E0 from atomic energies.
const emt_parameters *singleelement; ///< Element if only one element, 0 otherwise.
std::vector parameters; ///< The EMT parameters
const TinyDoubleMatrix *chi; ///< The Chi matrix of EMT.
int nelements; ///< The number of different elements in the simulation
/// Cutoff parameters (from EMTParameterProvider).
double rFermi, rNbCut;
/// Cutoff slope
double cutoffslope;
bool sigma2isvalid; ///< For consistency checks.
bool stresses_are_normalized; ///< For consistency checks.
bool initialized; ///< SetAtoms has been called.
/// Temporary data for the atoms
/// Each atom has two sigmas for each possible neighboring element.
/// So these are arrays (nelements long).
vector > sigma1;
vector > sigma2;
//@{
/// Each atom has a single Ec, Eas, radius etc.
vector Ec;
vector Eas;
vector Epot;
vector radius;
vector dEds;
vector force;
/// The stresses.
vector stress;
//@}
/// Temporary data
vector tmp_double;
vector ex2; // Data stored between calls to CalculateEnergiesAfterSigma (with and without actual energy calculation)
int nAtomsRes, nSizeRes;
///< If there are ghostatoms, some extra space is reserved for the arrays
/// The atomic numbers are translated into IDs, integers in [0, nelements-1]
vector id;
/* The buffers for batch processing */
// static const int BUFLEN; // The batch--buffer size. Avoid powers of 2
// static const int NMAXELEMENTS;
double unnormalizedstress[6]; ///< The total stress before division by volume
double unitcellvolume; ///< Unit cell volume stored by GetStresses
/// A structure of counters to check if recalculations are necessary.
struct {
int ids;
int nblist;
int sigma1;
int sigma2;
int beforeforces;
int energies;
int forces;
int stresses;
int fullstresses;
int fullstresses_mom;
} counters;
/// A structure of bools, to remember if recalculations are necessary.
struct {
bool ids;
bool nblist;
bool sigma1;
bool sigma2;
bool beforeforces;
bool energies;
bool forces;
bool stresses;
bool fullstresses;
bool fullstresses_mom;
} recalc;
/// A flag indicating that Atoms has been opened in a child class.
bool skip_begin;
};
} // end namespace
#endif // ! _EMT_H