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 <string.h>
  13 #include <sstream>
  14 #include <stdio.h>
  15 #include <stdarg.h>
  16
  17 // #define BZ_DEBUG
  18
  19 #include <blitz/array.h>        // http://www.oonumerics.org/blitz/
  20 #include <random/uniform.h>     // blitz Random Number Generator - Mersenne Twister type
  21 #include "SimpleIni.h"          // Ini File Parser
  22
  23 using namespace std;
  24 using namespace blitz;
  25 using namespace ranlib;
  26
  27
  28
  29
  30 enum Occupation {
  31         empty = 0,
  32         CO
  33 };
  34
  35 enum Atom {
  36         no_atom = 0,
  37         Pd = 2
  38 };
  39
  40 class Lattice {
  41 public:
  42         Array<Occupation, 2> sites;
  43         Array<Atom, 2> surface;
  44         int Latitude, Longitude;
  45         int UnitCellSizeX, UnitCellSizeY;
  46 };
  47
  48 class Directions {
  49 public:
  50         int x;
  51         int y;
  52 };
  53
  54 class Interaction_Parameter {
  55 public:
  56         string name;
  57         Array<Directions, 1> direction;
  58 };
  59
  60 void usage(const char *prog)
  61 {
  62         fprintf(stderr, "usage: %s structure-ini-file interactions-ini-file\n", prog);
  63         exit(EXIT_FAILURE);
  64 }
  65
  66 void panic(char *msg, ...)
  67 {
  68         va_list vargs;
  69
  70         va_start(vargs, msg);
  71         fprintf(stderr, "FATAL: ");
  72         vfprintf(stderr, msg, vargs);
  73         fprintf(stderr, "\n");
  74         exit(EXIT_FAILURE);
  75 }
  76
  77 const char *str_SI_Error(SI_Error rc)
  78 {
  79         switch (rc) {
  80         case SI_OK:             return "no error";
  81         case SI_UPDATED:        return "existing value updated";
  82         case SI_INSERTED:       return "new value inserted";
  83         case SI_FAIL:           return "generic failure";
  84         case SI_NOMEM:          return "out of memory";
  85         case SI_FILE:           return "file error";
  86         default:                return "(invalid error)";
  87         }
  88 };
  89
  90 void read_ini_structure(const char *input_file_name, Lattice &lattice)
  91 {
  92         CSimpleIni ini(false, true, false);
  93         const char *tmp;
  94
  95         SI_Error rc = ini.LoadFile(input_file_name);
  96         if (rc < 0)
  97                 panic("Failed loading ini file: %s", str_SI_Error(rc));
  98
  99         // Parse UnitCell
 100         const CSimpleIniA::TKeyVal * section = ini.GetSection("UnitCell");
 101         if (!section)
 102                 panic("Ini file is missing 'UnitCell' section.");
 103
 104 #define ASSIGN_KEY(section, key_name, key_member) \
 105         tmp = ini.GetValue(section, key_name); \
 106         if (!tmp) \
 107                 panic("%s definition is missing '%s'", section, key_name); \
 108         lattice.key_member = atoi(tmp)
 109
 110         ASSIGN_KEY("UnitCell", "SizeX", UnitCellSizeX);
 111         ASSIGN_KEY("UnitCell", "SizeY", UnitCellSizeY);
 112         ASSIGN_KEY("UnitCell", "Longitude", Longitude);
 113         ASSIGN_KEY("UnitCell", "Latitude", Latitude);
 114 #undef ASSIGN_KEY
 115
 116         if (lattice.UnitCellSizeX < 1 || lattice.UnitCellSizeY < 1 ||
 117                         lattice.Longitude < 1 || lattice.Latitude < 1)
 118                 panic("Unit cell has unreasonable size.");
 119
 120         // Building adsorbate sites array
 121         lattice.sites.resize(lattice.Longitude * lattice.UnitCellSizeX + 1,
 122                         lattice.Latitude * lattice.UnitCellSizeY + 1);
 123         lattice.sites = empty;
 124
 125         // Building surface sites array
 126         lattice.surface.resize(lattice.Longitude * lattice.UnitCellSizeX + 1,
 127                         lattice.Latitude * lattice.UnitCellSizeY + 1);
 128         lattice.surface = no_atom;
 129
 130         // Parse Adsorbates
 131         section = ini.GetSection("Surface");
 132         if (!section)
 133                 panic("Ini file is missing 'Surface' section. Maybe meaningful??");
 134
 135         CSimpleIniA::TNamesDepend Pd_positions;
 136         ini.GetAllValues("Surface", "Pd", Pd_positions);
 137         for (CSimpleIniA::TNamesDepend::const_iterator i = Pd_positions.begin();
 138                         i != Pd_positions.end(); ++i) {
 139                 double x_frac, y_frac;
 140                 if (sscanf(i->pItem, "%lf/%lf", &x_frac, &y_frac) != 2 ||
 141                                 x_frac < 0 || x_frac >= 1 ||
 142                                 y_frac < 0 || y_frac >= 1)
 143                         panic("Skrewed coordinates in '%s'", i->pItem);
 144                 lattice.surface((int)(lattice.UnitCellSizeX * x_frac),
 145                                 (int)(lattice.UnitCellSizeY * y_frac)) = Pd;
 146         }
 147
 148         // Replicate surface unit cell to full dimension
 149         for (int x = 0; x < lattice.Longitude * lattice.UnitCellSizeX + 1; ++x)
 150                 for (int y = 0; y < lattice.Latitude * lattice.UnitCellSizeY + 1; ++y)
 151                         lattice.surface(x, y) =
 152                                 lattice.surface((x % lattice.UnitCellSizeX),
 153                                                 (y % lattice.UnitCellSizeY));
 154
 155         // Parse Adsorbates
 156         section = ini.GetSection("Adsorbates");
 157         if (!section)
 158                 panic("Ini file is missing 'Adsorbates' section. Maybe meaningful??");
 159
 160         CSimpleIniA::TNamesDepend CO_positions;
 161         ini.GetAllValues("Adsorbates", "CO", CO_positions);
 162         for (CSimpleIniA::TNamesDepend::const_iterator i = CO_positions.begin();
 163                         i != CO_positions.end(); ++i) {
 164                 double x_frac, y_frac;
 165                 if (sscanf(i->pItem, "%lf/%lf", &x_frac, &y_frac) != 2 ||
 166                                 x_frac < 0 || x_frac >= 1 ||
 167                                 y_frac < 0 || y_frac >= 1)
 168                         panic("Skrewed coordinates in '%s'", i->pItem);
 169                 lattice.sites((int)(lattice.UnitCellSizeX * x_frac),
 170                                 (int)(lattice.UnitCellSizeY * y_frac)) = CO;
 171         }
 172
 173         // Replicate adsorbates unit cell to full dimension
 174         for (int x = 0; x < lattice.Longitude * lattice.UnitCellSizeX + 1; ++x)
 175                 for (int y = 0; y < lattice.Latitude * lattice.UnitCellSizeY + 1; ++y)
 176                         lattice.sites(x, y) =
 177                                 lattice.sites((x % lattice.UnitCellSizeX),
 178                                                 (y % lattice.UnitCellSizeY));
 179
 180 }
 181
 182
 183 void read_ini_interactions(const char *input_file_name, Array<Interaction_Parameter, 1> &interaction, int &total_interactions)
 184 {
 185         CSimpleIni ini(false, true, false);
 186
 187
 188         SI_Error rc = ini.LoadFile(input_file_name);
 189         if (rc < 0)
 190                 panic("Failed loading ini file: %s", str_SI_Error(rc));
 191
 192         // Parse Interactions Section
 193         const CSimpleIniA::TKeyVal * section = ini.GetSection("Interactions");
 194         if (!section)
 195                 panic("Interactions.ini file is missing 'Interactions' section.");
 196
 197         // First we need to get all keys and provide an interaction array of the correct size
 198         //total_interactions = 0;
 199
 200         CSimpleIniA::TNamesDepend Interaction_types;
 201         ini.GetAllKeys("Interactions", Interaction_types);
 202
 203         // Counting all key value pairs for each key type
 204         for (CSimpleIniA::TNamesDepend::const_iterator i = Interaction_types.begin();
 205                         i != Interaction_types.end(); ++i) {
 206                                 CSimpleIniA::TNamesDepend Interaction_parameters;
 207                                 ini.GetAllValues("Interactions", i->pItem, Interaction_parameters);
 208                                 total_interactions += (int) (Interaction_parameters.size());
 209                         };
 210
 211         cout << " total interactions after counting " << total_interactions << endl;
 212
 213         interaction.resize(total_interactions);
 214
 215         // Counter to access interactions array
 216         int interaction_counter = 0;
 217
 218         // constructor for directions array
 219         for (int j = 0; j < total_interactions; j++){
 220                 interaction(j).direction(1);
 221         };
 222
 223         cout << " gone through interactions array " << endl;
 224
 225
 226         for (CSimpleIniA::TNamesDepend::const_iterator i = Interaction_types.begin();
 227                         i != Interaction_types.end(); ++i) {
 228                                 // Then we need to get all values for this key
 229                                 CSimpleIniA::TNamesDepend Interaction_parameters;
 230                                 ini.GetAllValues("Interactions", i->pItem, Interaction_parameters);
 231                                 cout << " getting values for " << i->pItem << endl;
 232                         for (CSimpleIniA::TNamesDepend::const_iterator k = Interaction_parameters.begin();
 233                         k != Interaction_parameters.end(); ++k) {
 234
 235                         cout << "getting first token of " << k->pItem << endl;
 236
 237                                 int resize_factor = 1; // at least one possible interaction vector for each interaction parameter
 238                                 int x_vec;
 239                                 int y_vec;
 240
 241                                 char * p_key_value = NULL;
 242                                 char * pch = NULL;
 243
 244                                 p_key_value = (strdup(k->pItem));
 245                                 cout << p_key_value << endl;
 246                                 // char *p_key_value = &key_value; // pointer to working copy
 247                                 // cout << p_key_value << *p_key_value << endl << endl;
 248                                                                                 // pointer for string tokens
 249
 250                                 pch = strtok(p_key_value," :"); // get first string token
 251
 252                                 cout << "first token is: " << pch << endl << endl;
 253
 254                                 pch = strtok(NULL,"  :");
 255
 256                                 cout << "second token is: " << pch << endl;
 257
 258                                 while (pch != NULL) {                   // until not reached the end of the interaction parameter line
 259                                         if (sscanf(pch, "%i/%i", &x_vec, &y_vec) != 2) {        // check if it is a x/y pair
 260                                                 cout << "assigning name " << endl;
 261                                                 interaction(interaction_counter).name = *pch;           // no x/y pair, so it's name of interaction
 262                                                 pch = strtok(NULL,"  :");
 263                                         }
 264                                         else {          // it is an x/y pair
 265                                                 cout << "assigning direction " << endl;
 266                                                 interaction(interaction_counter).direction.resize(resize_factor); // resize directions array, provide memory
 267                                                 interaction(interaction_counter).direction((resize_factor-1)).x = x_vec;        // assign values
 268                                                 interaction(interaction_counter).direction((resize_factor-1)).y = y_vec;
 269                                                 resize_factor++; // increase directions counter for possible next interaction
 270                                                 pch = strtok(NULL," :"); // get next token
 271                                                 };
 272                                         interaction_counter++;          // reached end of interaction parameter line, increase interactions counter
 273                                         resize_factor = 1;                      // reset directions counter for next interaction parameter
 274                                 };
 275
 276                         }; // end for loop for each interaction parameter
 277                         }; // end for loop for each interaction type
 278
 279 // some output for checking interaction array content
 280         for (int i=0; i < total_interactions; i++){
 281                 cout << "i= " << i << endl;
 282                 cout << interaction(i).name << endl;
 283                 cout << interaction(i).direction(0).x << endl;
 284         };
 285
 286  cout << "leaving interaction reading !! " << endl;
 287  cout << "total interactions " << total_interactions << endl;
 288 }
 289
 290
 291 int main(int argc, char* argv[])
 292 {
 293         Lattice lattice;
 294         Array<Interaction_Parameter, 1> interaction;
 295         int total_interactions;
 296         total_interactions = 0;
 297
 298
 299
 300         if (argc != 3)
 301                 usage(*argv);
 302         read_ini_structure(argv[1], lattice);
 303
 304         cout << "lattice:" << endl << lattice.sites;
 305         cout << "underlying surface:" << endl << lattice.surface;
 306
 307
 308         cout << "total interactions before interaction ini reading " << total_interactions << endl;
 309
 310         read_ini_interactions(argv[2], interaction, total_interactions);
 311
 312         cout << " now content output " << endl;
 313         cout << "total interactions  " << total_interactions << endl;
 314         /*
 315         // some output for checking interaction array content
 316         for (int i=0; i < 5; i++){
 317                 cout << "i= " << i << endl;
 318                 // cout << interaction(i).name << endl;
 319                 cout << interaction(i).direction(0).x << endl;
 320         };
 321         */
 322
 323         exit(EXIT_SUCCESS);
 324 }