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