?  search:
summary | log | graphiclog1 | graphiclog2 | commit | commitdiff | tree | refs | edit | fork
blame | history | raw | HEAD
added control output of variable total_interactions before reading interactions.ini
[cluster_expansion.git] / main.cpp
blob34bfc1546a605b47a3051019440b2e20b8bf9225
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 string key_name;
232 key_name = i->pItem;
233 // Checking what kind of interaction it is, choosing correct reading routine
234 if (key_name == "pair" ){
235 int x_vec;
236 int y_vec;
237 char Interaction_name;
238 cout << "interaction counter before " << interaction_counter << endl;
239 if (sscanf(k->pItem, "%i/%i \"%s\"", &x_vec, &y_vec, &Interaction_name) != 3 )
240 panic("Skrewed Pair Interaction description '%s'", k->pItem);
241 interaction(interaction_counter).name = Interaction_name;
242 interaction(interaction_counter).direction.resize(1);
243 interaction(interaction_counter).direction(0).x = x_vec;
244 interaction(interaction_counter).direction(0).y = y_vec;
245 interaction_counter++;
246 cout << "interaction counter after pairs " << interaction_counter << endl;
247 }
248 else if (key_name == "trio"){
249 int x_vec;
250 int x1_vec;
251 int y_vec;
252 int y1_vec;
253 char Interaction_name;
254 cout << "interaction counter before trio" << interaction_counter << endl;
255 if (sscanf(k->pItem, "%i/%i:%i/%i \"%s\"", &x_vec, &y_vec, &x1_vec, &y1_vec, &Interaction_name) != 5 )
256 panic("Skrewed Trio Interaction description '%s'", k->pItem);
257 interaction(interaction_counter).name = Interaction_name;
258 interaction(interaction_counter).direction.resize(2);
259 interaction(interaction_counter).direction(0).x = x_vec;
260 interaction(interaction_counter).direction(0).y = y_vec;
261 interaction(interaction_counter).direction(1).x = x1_vec;
262 interaction(interaction_counter).direction(1).y = y1_vec;
263 interaction_counter++;
264 cout << "interaction counter after trios " << interaction_counter << endl;
265 }
266 else if (key_name == "quattro"){
267 int x_vec;
268 int x1_vec;
269 int x2_vec;
270 int y_vec;
271 int y1_vec;
272 int y2_vec;
273 char Interaction_name;
274 if (sscanf(k->pItem, "%i/%i:%i/%i:%i/%i \"%s\"", &x_vec, &y_vec, &x1_vec, &y1_vec,
275 &x2_vec, &y2_vec, &Interaction_name) != 7 )
276 panic("Skrewed Quattro Interaction description '%s'", k->pItem);
277 interaction(interaction_counter).name = Interaction_name;
278 interaction(interaction_counter).direction.resize(3);
279 interaction(interaction_counter).direction(0).x = x_vec;
280 interaction(interaction_counter).direction(0).y = y_vec;
281 interaction(interaction_counter).direction(1).x = x1_vec;
282 interaction(interaction_counter).direction(1).y = y1_vec;
283 interaction(interaction_counter).direction(2).x = x2_vec;
284 interaction(interaction_counter).direction(2).y = y2_vec;
285 interaction_counter++;
286 }
287 else if (key_name == "quinto"){
288 int x_vec;
289 int x1_vec;
290 int x2_vec;
291 int x3_vec;
292 int y_vec;
293 int y1_vec;
294 int y2_vec;
295 int y3_vec;
296 char Interaction_name;
297 if (sscanf(k->pItem, "%i/%i:%i/%i:%i/%i:%i/%i \"%s\"", &x_vec, &y_vec, &x1_vec, &y1_vec,
298 &x2_vec, &y2_vec, &x3_vec, &y3_vec, &Interaction_name) != 9 )
299 panic("Skrewed Quinto Interaction description '%s'", k->pItem);
300 interaction(interaction_counter).name = Interaction_name;
301 interaction(interaction_counter).direction.resize(4);
302 interaction(interaction_counter).direction(0).x = x_vec;
303 interaction(interaction_counter).direction(0).y = y_vec;
304 interaction(interaction_counter).direction(1).x = x1_vec;
305 interaction(interaction_counter).direction(1).y = y1_vec;
306 interaction(interaction_counter).direction(2).x = x2_vec;
307 interaction(interaction_counter).direction(2).y = y2_vec;
308 interaction(interaction_counter).direction(3).x = x3_vec;
309 interaction(interaction_counter).direction(3).y = y3_vec;
310 interaction_counter++;
311 };
312 };
313 };
314 cout << "leaving interaction reading !! " << endl;
315 cout << "total interactions " << total_interactions << endl;
316 }
317
318
319 int main(int argc, char* argv[])
320 {
321 Lattice lattice;
322 Array<Interaction_Parameter, 1> interaction;
323 int total_interactions;
324 total_interactions = 0;
325
326
327
328 if (argc != 3)
329 usage(*argv);
330 read_ini_structure(argv[1], lattice);
331
332 cout << "lattice:" << endl << lattice.sites;
333 cout << "underlying surface:" << endl << lattice.surface;
334
335
336 cout << "total interactions before interaction ini reading " << total_interactions << endl;
337
338 read_ini_interactions(argv[2], interaction, total_interactions);
339
340 cout << " now content output " << endl;
341 cout << " total interactions are " << total_interactions;
342
343 // some output for checking interaction array content
344 for (int i=0; i < total_interactions; i++){
345 cout << "i= " << i << endl;
346 cout << interaction(i).name << endl;
347 cout << interaction(i).direction(0).x << endl;
348 };
349
350
351 exit(EXIT_SUCCESS);
352 }
Cluster Expansion for CO lateral interactions on Pd(100)