main.cpp

   1 /*
   2  * Determination of lateral interaction parameters on a c(2x2) square lattice
   3  * using linear combinations to express each interaction and then counting
   4  * them in terms of cluster-expansion coefficients
   5  *
   6  * Michael Rieger, FHI, 2008
   7  *
   8  */
   9
  10 #include <iostream>
  11 #include <string>
  12 #include <sstream>
  13 #include <stdio.h>
  14 #include <stdarg.h>
  15
  16 #include <blitz/array.h>        // http://www.oonumerics.org/blitz/
  17 #include <random/uniform.h>     // blitz Random Number Generator - Mersenne Twister type
  18 #include "SimpleIni.h"          // Ini File Parser
  19
  20 using namespace std;
  21 using namespace blitz;
  22 using namespace ranlib;
  23
  24
  25 enum Occupation {
  26         empty = 0,
  27         CO
  28 };
  29
  30 enum Atom {
  31         no_atom = 0,
  32         Pd = 2
  33 };
  34
  35 enum Pair_Interaction_Parameters {
  36         no_pair_interaction = 0,
  37         Vp1,
  38         Vp2,
  39         Vp2across,
  40         Vp3,
  41         Vp4,
  42         Vp4across,
  43         Vp5,
  44         Vp5across,
  45         Vp6,
  46         max_pair_interactions
  47 };
  48
  49 enum Trio_Interaction_Parameters {
  50         no_trio_interaction = 0,
  51         Vtr,
  52         Vtr112,
  53         Vtr112across,
  54         Vtr2,
  55         Vtr2across,
  56         Vtr22across3,
  57         Vtr214,
  58         Vtr214across,
  59         Vtr3,
  60         Vtr13,
  61         Vtr335,
  62         Vtr335across,
  63         Vtr223across,
  64         Vtr134,
  65         Vtr237,
  66         Vtr237across,
  67         max_trio_interactions
  68 };
  69
  70 class Lattice {
  71 public:
  72         Array<Occupation, 2> sites;
  73         Array<Atom, 2> surface;
  74         int Latitude, Longitude;
  75         int UnitCellSizeX, UnitCellSizeY;
  76 };
  77
  78 struct directions {
  79         int x;
  80         int y;
  81 };
  82
  83 // direction.x = etc..
  84
  85 struct interaction_definition {
  86         string name;
  87         Array<directions, 1> direction;
  88 };
  89
  90
  91 // interaction_parameter.name and interaction_parameter.interaction_vector(direction.x,1)
  92 // and interaction_parameter.interaction_vector(direction.y,1)
  93
  94 Array<interaction_definition, 1> Interactions;
  95
  96
  97 void usage(const char *prog)
  98 {
  99         fprintf(stderr, "usage: %s structure-ini-file\n", prog);
 100         exit(EXIT_FAILURE);
 101 }
 102
 103 void panic(char *msg, ...)
 104 {
 105         va_list vargs;
 106
 107         va_start(vargs, msg);
 108         fprintf(stderr, "FATAL: ");
 109         vfprintf(stderr, msg, vargs);
 110         fprintf(stderr, "\n");
 111         exit(EXIT_FAILURE);
 112 }
 113
 114 const char *str_SI_Error(SI_Error rc)
 115 {
 116         switch (rc) {
 117         case SI_OK:             return "no error";
 118         case SI_UPDATED:        return "existing value updated";
 119         case SI_INSERTED:       return "new value inserted";
 120         case SI_FAIL:           return "generic failure";
 121         case SI_NOMEM:          return "out of memory";
 122         case SI_FILE:           return "file error";
 123         default:                return "(invalid error)";
 124         }
 125 };
 126
 127 void read_ini_structure(const char *input_file_name, Lattice &lattice)
 128 {
 129         CSimpleIni ini(false, true, false);
 130         const char *tmp;
 131
 132         SI_Error rc = ini.LoadFile(input_file_name);
 133         if (rc < 0)
 134                 panic("Failed loading ini file: %s", str_SI_Error(rc));
 135
 136         // Parse UnitCell
 137         const CSimpleIniA::TKeyVal * section = ini.GetSection("UnitCell");
 138         if (!section)
 139                 panic("Ini file is missing 'UnitCell' section.");
 140
 141 #define ASSIGN_KEY(section, key_name, key_member) \
 142         tmp = ini.GetValue(section, key_name); \
 143         if (!tmp) \
 144                 panic("%s definition is missing '%s'", section, key_name); \
 145         lattice.key_member = atoi(tmp)
 146
 147         ASSIGN_KEY("UnitCell", "SizeX", UnitCellSizeX);
 148         ASSIGN_KEY("UnitCell", "SizeY", UnitCellSizeY);
 149         ASSIGN_KEY("UnitCell", "Longitude", Longitude);
 150         ASSIGN_KEY("UnitCell", "Latitude", Latitude);
 151 #undef ASSIGN_KEY
 152
 153         if (lattice.UnitCellSizeX < 1 || lattice.UnitCellSizeY < 1 ||
 154                         lattice.Longitude < 1 || lattice.Latitude < 1)
 155                 panic("Unit cell has unreasonable size.");
 156
 157         // Building adsorbate sites array
 158         lattice.sites.resize(lattice.Longitude * lattice.UnitCellSizeX + 1,
 159                         lattice.Latitude * lattice.UnitCellSizeY + 1);
 160         lattice.sites = empty;
 161
 162         // Building surface sites array
 163         lattice.surface.resize(lattice.Longitude * lattice.UnitCellSizeX + 1,
 164                         lattice.Latitude * lattice.UnitCellSizeY + 1);
 165         lattice.surface = no_atom;
 166
 167         // Parse Adsorbates
 168         section = ini.GetSection("Surface");
 169         if (!section)
 170                 panic("Ini file is missing 'Surface' section. Maybe meaningful??");
 171
 172         CSimpleIniA::TNamesDepend Pd_positions;
 173         ini.GetAllValues("Surface", "Pd", Pd_positions);
 174         for (CSimpleIniA::TNamesDepend::const_iterator i = Pd_positions.begin();
 175                         i != Pd_positions.end(); ++i) {
 176                 double x_frac, y_frac;
 177                 if (sscanf(i->pItem, "%lf/%lf", &x_frac, &y_frac) != 2 ||
 178                                 x_frac < 0 || x_frac >= 1 ||
 179                                 y_frac < 0 || y_frac >= 1)
 180                         panic("Skrewed coordinates in '%s'", i->pItem);
 181                 lattice.surface((int)(lattice.UnitCellSizeX * x_frac),
 182                                 (int)(lattice.UnitCellSizeY * y_frac)) = Pd;
 183         }
 184
 185         // Replicate surface unit cell to full dimension
 186         for (int x = 0; x < lattice.Longitude * lattice.UnitCellSizeX + 1; ++x)
 187                 for (int y = 0; y < lattice.Latitude * lattice.UnitCellSizeY + 1; ++y)
 188                         lattice.surface(x, y) =
 189                                 lattice.surface((x % lattice.UnitCellSizeX),
 190                                                 (y % lattice.UnitCellSizeY));
 191
 192         // Parse Adsorbates
 193         section = ini.GetSection("Adsorbates");
 194         if (!section)
 195                 panic("Ini file is missing 'Adsorbates' section. Maybe meaningful??");
 196
 197         CSimpleIniA::TNamesDepend CO_positions;
 198         ini.GetAllValues("Adsorbates", "CO", CO_positions);
 199         for (CSimpleIniA::TNamesDepend::const_iterator i = CO_positions.begin();
 200                         i != CO_positions.end(); ++i) {
 201                 double x_frac, y_frac;
 202                 if (sscanf(i->pItem, "%lf/%lf", &x_frac, &y_frac) != 2 ||
 203                                 x_frac < 0 || x_frac >= 1 ||
 204                                 y_frac < 0 || y_frac >= 1)
 205                         panic("Skrewed coordinates in '%s'", i->pItem);
 206                 lattice.sites((int)(lattice.UnitCellSizeX * x_frac),
 207                                 (int)(lattice.UnitCellSizeY * y_frac)) = CO;
 208         }
 209
 210         // Replicate adsorbates unit cell to full dimension
 211         for (int x = 0; x < lattice.Longitude * lattice.UnitCellSizeX + 1; ++x)
 212                 for (int y = 0; y < lattice.Latitude * lattice.UnitCellSizeY + 1; ++y)
 213                         lattice.sites(x, y) =
 214                                 lattice.sites((x % lattice.UnitCellSizeX),
 215                                                 (y % lattice.UnitCellSizeY));
 216
 217 }
 218
 219 int main(int argc, char* argv[])
 220 {
 221         Lattice lattice;
 222         if (argc != 2)
 223                 usage(*argv);
 224         read_ini_structure(argv[1], lattice);
 225
 226         cout << "lattice:" << endl << lattice.sites;
 227         cout << "underlying surface:" << endl << lattice.surface;
 228
 229         exit(EXIT_SUCCESS);
 230 }