Structure.cpp

   1 #include <stdexcept>
   2
   3 #include <blitz/array.h>
   4 #include <random/uniform.h>                             // blitz RNG Mersenne Twister
   5 #include <set>
   6 using namespace std;
   7 using namespace blitz;
   8 using namespace ranlib;
   9
  10 #include "Lattice.h"
  11 #include "Interactions.h"
  12 #include "Structure.h"
  13
  14 /*
  15 #define WANT_STREAM                             // include.h will get stream fns
  16 //#define WANT_MATH                             // include.h will get math fns
  17 #define USING_DOUBLE                                            // defines newmat with double precision
  18                                                                                         // - we can maybe use this tyedefinition of Real type in newmat to get even more precision with libgmp.a
  19 // newmatap.h will get include.h
  20
  21 #include "newmat/newmatap.h"                // need matrix applications
  22 #include "newmat/newmatio.h"                // need matrix output routines
  23 #ifdef use_namespace
  24 using namespace NEWMAT; // access NEWMAT namespace
  25 #endif
  26 */
  27
  28
  29
  30 void Structures::printEnergies(Structures& structures) {
  31         int counter = 1;
  32         double meanSquareDifference = 0.0;
  33
  34         for (Structures::iterator structure = structures.begin(); structure != structures.end(); structure++){
  35                 double deltaE = ((*structure)->dftEnergy - (*structure)->onSiteEnergy);
  36                 cout << "Structure Nr: " << counter << " got: \n" <<
  37                 "\t CE  energy of: " << (*structure)->getEnergy() << "\n" <<
  38                 "\t - DFT energy of: " << (*structure)->dftEnergy << "\n" <<
  39                 "\t onSite energy: " << (*structure)->onSiteEnergy << "\n" <<
  40                 "\t DFT - onSite: " << deltaE << "\n" <<
  41                 "\t delta(CE - DFT): " <<  ((*structure)->getEnergy() - deltaE)   << "\n";
  42                 meanSquareDifference += ((*structure)->getEnergy() - deltaE) * ((*structure)->getEnergy() - deltaE);
  43                 cout << " meanSquareDifference: " << meanSquareDifference << "\n";
  44                 counter++;
  45         };
  46
  47         double CVScore = sqrt( meanSquareDifference / ((double) counter));
  48
  49         cout << " finally the CV Score for this: " << CVScore << "\n";
  50
  51         cout << " now showing the power of combinations.hpp: \n";
  52
  53         Structures::iterator referenceStructure = structures.begin();
  54
  55         // just a test to see, if i can call the combinations routine here
  56         //(*referenceStructure)->interactions.generateCombinations();
  57
  58         return;
  59 }
  60
  61
  62 double Structure::getEnergy() {
  63         // FIXME: check adsorbate species and get correct on-site energy contribution
  64
  65         // get interaction multiplicities for current configuration
  66         // lattice->assessInteractions(interactions);
  67
  68         double energy = 0.0;
  69         // energy += onSiteEnergy; // TODO : adsorbate_on_site_energy(species and site specific!)
  70                                                         // so far only sum value for all species present solved in input file
  71         /* no need for this yet, since onSiteEnergy gives sum value
  72         for (LatticeLayer::iterator i = lattice->adsorbates.begin(); i
  73                         != lattice->adsorbates.end(); i++) {
  74                 if (lattice->adsorbates(i.position()) == empty)
  75                         continue;
  76                 // implement species checking CO or O or whatever
  77                 };
  78         */
  79         for (Interactions::iterator interaction = interactions.begin(); interaction
  80                         != interactions.end(); interaction++) {
  81                 // cout << "calculating energy contribution of " << endl
  82                 //      << interaction->name << " with " << interaction->multiplicity << "*" << interaction->energy << endl;
  83                 int interactionMultiplier = 0; // to take into account that we fitted all interaction energies normalized to 1 atom
  84                 switch (interaction->directions.size()) {
  85                         case (1) : { interactionMultiplier = 2; break; }; // pair interaction
  86                         case (2) : { interactionMultiplier = 3; break; }; // trio interaction
  87                         case (3) : { interactionMultiplier = 4; break; }; // quattro interaction
  88                         case (4) : { interactionMultiplier = 5; break; }; // quinto interaction
  89                 };
  90                 energy += (interaction->energy * interaction->multiplicity); // * interactionMultiplier);
  91         };
  92
  93         //energy *= lattice->nrCO;
  94         energy += onSiteEnergy;
  95                         // cout << "total structure energy is: " << energy << endl;
  96         //energy = energy * (-1); // changing sign of energy - TODO: why is that necessary?
  97         return energy;
  98 }
  99
 100