EMT.h

Go to the documentation of this file.

 
#ifndef _EMT_H
#define _EMT_H
 
#include "EMTParameterProvider.h"
#include "Potential.h"
#include "TinyMatrix.h"
#include <vector>
using std::vector;
 
class AtomsBase;
class Atoms;
class GhostAtoms;
class NeighborList;
struct emt_parameters;


class EMT : public AsapPotential {
public:
  EMT(EMTParameterProvider *prov = 0);

  virtual ~EMT();
  void SetSubtractE0(bool subtractE0) { this->subtractE0 = subtractE0; }
 
  virtual void SetAtoms(Atoms *atoms);
  virtual const Vec *GetCartesianForces();
  virtual const double *GetPotentialEnergies();
  virtual double GetPotentialEnergy();
  virtual const symTensor *GetStresses(const Vec *momenta = 0);
  virtual void GetStress(double stress[6], const Vec *momenta = 0);
  virtual void CheckNeighborLists();
 
  virtual double GetCutoffRadius() const { return rNbCut; }
  virtual double GetLatticeConstant() const;
  virtual int GetNumberOfAtoms() const { return nAtoms; }

  virtual void SetContinuumElement(int z);
  virtual void UpdateSuperCell(const SuperCell *newSuperCell);
 
  // Quasicontinuum stuff
  virtual double CalculateLatticeEnergy(const Vec a[3]);
  virtual void CalculateDerivatives(const Vec a[3], double dEdaDota[6]);
  virtual double GetData() const { return dataSigma1; }

  const double *GetSigma1(int n) { return &sigma1[n][0]; }
  const double *GetSigma2(int n) { return &sigma2[n][0]; }

  void PrintParameters();

  NeighborList *GetNeighborList() const { return nblist; }
 
protected:
  virtual void InitParameters();
  virtual void Allocate();
  virtual void AllocateStress();
  virtual void CalculateIDs();
  virtual void CalculateSigmas(int calculatesigma2 = 1);

  virtual void CalculateEnergiesAfterSigmas(double *Epot = 0);

  virtual void CalculateForcesAfterEnergies(Vec *forces,
                                            double stress[][6] = 0);

  virtual void CalculateForcesAfterEnergiesSingle(Vec *forces,
                                                  double stress[][6] = 0);

  virtual void CalculateStuff(const Vec a[3], double &sigma1, double &sigma2,
                              double *diffs2, double *wght, double *dsigma1,
                              double *dsigma2);
 
private: // Methods.

  void sigma_batch(int *self, int *other, Vec rnb[], double *sq_dist, int zs,
                   int zo, int n, int calculatesigma2);


  void force_batch(int *self, int *other, Vec rnb[], double *sq_dist,
                   double dEdss[], double dEdso[], int zs, int zo, int n,
                   Vec *force, double (*stress)[6]);
 
protected:                        // Data
  Atoms *atoms;                   
  GhostAtoms *ghostatoms;         
  int nAtoms;                     
  int nSize;                      
  NeighborList *nblist;           
  EMTParameterProvider *provider; 
  bool ownProvider;               
 
  bool subtractE0; 
 
  std::vector<const emt_parameters *> parameters; 
  const emt_parameters *continuumelement;
  const TinyDoubleMatrix *chi; 
  int nelements; 
  double rFermi, rNbCut;
  double cutoffslope;
  //  int maxleaflen;         // Lenght of the longest leaf.
  int sigma2isvalid; 


  vector<vector<double>> sigma1;
  vector<vector<double>> sigma2;


  vector<double> Ec;
  vector<double> Eas;
  vector<double> potentialenergy;
  vector<double> radius;
  vector<double> dEds;
  vector<Vec> force;

  vector<double> stress;
 
  int nAtomsRes, nSizeRes;

  vector<int> 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]; 
  double totalvolume;           
 
  int nHalfNeighbors; 
  int (*coef)[6];     
  double dataSigma1;  

  struct {
    int ids;
    int nblist;
    int sigma1;
    int sigma2;
    int beforeforces;
    int energies;
    int forces;
    int stresses;
    int fullstresses;
  } counters;
};
 
#endif // ! _EMT_H