1 // EMT.h -- Implements the EMT potential. 2 // -*- c++ -*- 3 // 4 // Copyright (C) 2001-2011 Jakob Schiotz and Center for Individual 5 // Nanoparticle Functionality, Department of Physics, Technical 6 // University of Denmark. Email: schiotz@fysik.dtu.dk 7 // 8 // This file is part of Asap version 3. 9 // Asap is released under the GNU Lesser Public License (LGPL) version 3. 10 // However, the parts of Asap distributed within the OpenKIM project 11 // (including this file) are also released under the Common Development 12 // and Distribution License (CDDL) version 1.0. 13 // 14 // This program is free software: you can redistribute it and/or 15 // modify it under the terms of the GNU Lesser General Public License 16 // version 3 as published by the Free Software Foundation. Permission 17 // to use other versions of the GNU Lesser General Public License may 18 // granted by Jakob Schiotz or the head of department of the 19 // Department of Physics, Technical University of Denmark, as 20 // described in section 14 of the GNU General Public License. 21 // 22 // This program is distributed in the hope that it will be useful, 23 // but WITHOUT ANY WARRANTY; without even the implied warranty of 24 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 25 // GNU General Public License for more details. 26 // 27 // You should have received a copy of the GNU General Public License 28 // and the GNU Lesser Public License along with this program. If not, 29 // see <http://www.gnu.org/licenses/>. 30 31 #ifndef _EMT_H 32 #define _EMT_H 33 34 #include "AsapPython.h" 35 #include "Asap.h" 36 #include "Atoms.h" 37 #include "Vec.h" 38 #include "SymTensor.h" 39 #include "Potential.h" 40 #include "EMTParameterProvider.h" 41 #include "TinyMatrix.h" 42 #include <vector> 43 using std::vector; 44 45 namespace ASAPSPACE { 46 47 class NeighborLocator; 48 struct emt_parameters; 49 50 /// The Effective Medium Theory (EMT) potential. 51 52 /// The Effective Medium Theory (EMT) as originally implemented in Per 53 /// Stoltze's ARTwork program, and documented in 54 /// - K. W. Jacobsen, P. Stoltze and J. K. Norskov, 55 /// Surf. Sci. vol. 366, p. 394-402 (1996). 56 class EMT : public Potential 57 { 58 public: 59 /// Create an EMT potential optionally using a parameter provider. 60 EMT(PyObject *self, PyObject *prov, int verbose); 61 62 /// Delete the EMT potential. 63 virtual ~EMT(); 64 GetName()65 virtual string GetName() const {return "EMT";} 66 virtual string GetRepresentation() const; 67 68 /// Whether we subtract E0 from atomic energies. 69 /// 70 /// Defines the zero of potential energy; if we don't subtract, zero 71 /// corresponds to infinite seperation. If we do, it is the bulk fcc 72 /// crystal at the equilibrium lattice constant. SetSubtractE0(bool subtractE0)73 void SetSubtractE0(bool subtractE0) { this->subtractE0 = subtractE0; } 74 75 virtual void SetAtoms(PyObject *pyatoms, Atoms* accessobj = NULL); 76 virtual const vector<Vec> &GetForces(PyObject *pyatoms); 77 virtual const vector<double> &GetPotentialEnergies(PyObject *pyatoms); 78 virtual double GetPotentialEnergy(PyObject *pyatoms); 79 virtual const vector<SymTensor> &GetVirials(PyObject *pyatoms); 80 /// Check if neighbor lists need an update. 81 /// 82 /// Return true if the neigbor list needs an update. In parallel 83 /// simulations, the atoms will cause a communication so it returns 84 /// true if an update is needed on any processor. In parallel 85 /// simulations, it also instructs the atoms to communicate their 86 /// update status. 87 virtual bool CheckNeighborList(); 88 /// Update the neighbor lists 89 virtual void UpdateNeighborList(); 90 91 virtual bool CalcReq_Energy(PyObject *pyatoms); 92 virtual bool CalcReq_Forces(PyObject *pyatoms); 93 virtual bool CalcReq_Virial(PyObject *pyatoms); 94 virtual bool CalcReq_Virials(PyObject *pyatoms); 95 GetCutoffRadius()96 virtual double GetCutoffRadius() const {return rNbCut;} 97 virtual double GetLatticeConstant() const; GetNumberOfAtoms()98 virtual int GetNumberOfAtoms() const {return nAtoms;} 99 100 /// Return a pointer to the EMT "density" sigma1. GetSigma1()101 const vector<vector<double> > &GetSigma1() {return sigma1;} 102 103 /// Return a pointer to the EMT "density" sigma2. GetSigma2()104 const vector<vector<double> > &GetSigma2() {return sigma2;} 105 106 /// Print the EMT parameters. 107 void PrintParameters(); 108 109 /// Return the neighbor list. 110 111 /// Return a BORROWED reference to the Python object containing 112 /// neighbor list for the 113 /// potential, if this type of potential supports it, and if it is 114 /// defined now. Otherwise, return NULL without setting a Python 115 /// error. GetNeighborList()116 PyObject *GetNeighborList() const {return nblist_obj;} 117 118 /// This potential can be used in parallel simulations Parallelizable()119 virtual bool Parallelizable() const {return true;} 120 121 /// Mark that the atoms have changed the boundary conditions. 122 virtual void BoundaryConditionsChanged(); 123 124 /// Print memory usage 125 virtual long PrintMemory() const; 126 127 /// Return the atomic volumes that should be used to calculate the stresses 128 virtual void GetAtomicVolumes(vector<double> &volumes); 129 130 protected: 131 /// Initialization of the EMT parameters. 132 virtual void InitParameters(); 133 /// (Re)allocate storage for forces, energies and intermediate results. 134 virtual void Allocate(); 135 /// (Re)allocate storage for stresses 136 virtual void AllocateStress(); 137 /// Calculate type numbers from the atomic numbers. 138 virtual void CalculateIDs(); 139 140 /// Calculate sigma1 and perhaps sigma2. 141 virtual void CalculateSigmas(bool calculatesigma2); 142 143 /// Calculate energies 144 virtual void CalculateEnergies(); 145 146 /// Calculate energies once sigma1 is known 147 148 /// If Epot is NULL, energies are not calculated, only derivatives 149 /// needed for the force. 150 virtual void CalculateEnergiesAfterSigmas(bool calcEpot); 151 152 /// Calculate the forces 153 virtual void CalculateForces(); 154 155 /// Calculate forces in a multicomponent system. 156 virtual void CalculateForcesAfterEnergies(); 157 158 /// Calculate forces in a system with only one element. 159 virtual void CalculateForcesAfterEnergiesSingle(); 160 161 /// Calculate stresses. 162 virtual void CalculateVirials(); 163 164 /// Allocate the neighbor list. 165 virtual void CreateNeighborList(); 166 167 /// sigma_batch does the hard work in CalculateSigmas(). 168 virtual void sigma_batch(int *self, int *other, 169 double *sq_dist, int zs, int zo, int n, 170 bool calculatesigma2, bool partialupdate = false); 171 172 /// force_batch does the hard work in CalculateForcesAfterEnergy(). 173 virtual void force_batch(const int *self, const int *other, const Vec rnb[], 174 const double sq_dist[], const double dEdss[], 175 const double dEdso[], int zs, int zo, int n); 176 void distribute_force(const int *self, const int *other, 177 const double *df, const Vec *rnb, int n); 178 179 protected: // Data 180 bool ghostatoms; ///< True if number of atoms is likely to vary (atoms have ghosts). 181 int nAtoms; ///< The number of (real) atoms. 182 int nSize; ///< Number of atoms including ghost atoms. 183 NeighborLocator *nblist; ///< The neighborlist object. 184 PyObject *nblist_obj; 185 double driftfactor; ///< Drift factor for the neighbor list. 186 EMTParameterProvider *provider; ///< The source of the EMT parameters 187 PyObject *provider_obj; ///< The Python object behind provider. 188 bool subtractE0; ///< Whether we subtract E0 from atomic energies. 189 const emt_parameters *singleelement; ///< Element if only one element, 0 otherwise. 190 std::vector<const emt_parameters *> parameters; ///< The EMT parameters 191 const TinyDoubleMatrix *chi; ///< The Chi matrix of EMT. 192 int nelements; ///< The number of different elements in the simulation 193 /// Cutoff parameters (from EMTParameterProvider). 194 double rFermi, rNbCut; 195 /// Cutoff slope 196 double cutoffslope; 197 bool sigma2isvalid; ///< For consistency checks. 198 bool initialized; ///< SetAtoms has been called. 199 bool always_fullnblist; ///< Allways use full neighbor lists (beyond minimum-image). 200 201 /// Temporary data for the atoms 202 203 /// Each atom has two sigmas for each possible neighboring element. 204 /// So these are arrays (nelements long). 205 vector<vector<double> > sigma1; 206 vector<vector<double> > sigma2; 207 208 //@{ 209 /// Each atom has a single Ec, Eas, radius etc. 210 vector<double> Ec; 211 vector<double> Eas; 212 vector<double> Epot; 213 vector<double> radius; 214 vector<double> dEds; 215 vector<Vec> force; 216 /// The stresses. 217 vector<SymTensor> virials; 218 //@} 219 220 /// Temporary data 221 vector<double> tmp_double; 222 vector<double> ex2; // Data stored between calls to CalculateEnergiesAfterSigma (with and without actual energy calculation) 223 224 int nAtomsRes, nSizeRes; 225 ///< If there are ghostatoms, some extra space is reserved for the arrays 226 227 /// The atomic numbers are translated into IDs, integers in [0, nelements-1] 228 vector<int> id; 229 230 /* The buffers for batch processing */ 231 // static const int BUFLEN; // The batch--buffer size. Avoid powers of 2 232 // static const int NMAXELEMENTS; 233 234 /// A structure of counters to check if recalculations are necessary. 235 struct { 236 int ids; 237 int nblist; 238 int sigma1; 239 int sigma2; 240 int beforeforces; 241 int energies; 242 int forces; 243 int virials; 244 } counters; 245 246 /// A structure of bools, to remember if recalculations are necessary. 247 struct { 248 bool ids; 249 bool nblist; 250 bool sigma1; 251 bool sigma2; 252 bool beforeforces; 253 bool energies; 254 bool forces; 255 bool virials; 256 } recalc; 257 258 /// A flag indicating that Atoms has been opened in a child class. 259 bool skip_begin; 260 }; 261 262 } // end namespace 263 264 #endif // ! _EMT_H 265