Interactions.cpp

   1 #include <stdexcept>
   2 #include <string>
   3 #include <iostream>
   4
   5
   6 #include "Interactions.h"
   7 #include "SimpleIni.h"
   8 #include "Combination.hpp" // template class to generate all combinations of a vector - we need this for the interactions
   9
  10 using namespace std;
  11 using namespace boost; // needed for using the functions of combinations.hpp
  12 using namespace detail; // see above
  13
  14 ostream& operator<<(ostream& output, const Interactions& interactions) {
  15         for (Interactions::const_iterator interaction = interactions.begin(); interaction
  16                         != interactions.end(); interaction++) {
  17                 output << (interaction->name.empty() ? "(no name)" : interaction->name)
  18                                 << " [M:" << interaction->multiplicity << ", E:"
  19                                 << interaction->energy << "]\t\t: ";
  20                 for (Directions::const_iterator direction =
  21                                 interaction->directions.begin(); direction
  22                                 != interaction->directions.end(); direction++) {
  23                         output << "(" << direction->x << "/" << direction->y << ") ";
  24                 }
  25                 output << endl;
  26                 for (vector<Directions>::const_iterator directions =
  27                                 interaction->symmetryDirections.begin(); directions
  28                                 != interaction->symmetryDirections.end(); directions++) {
  29                         output << "\t\t\t";
  30                         for (Directions::const_iterator direction = directions->begin(); direction
  31                                         != directions->end(); direction++) {
  32                                 output << "(" << direction->x << "/" << direction->y
  33                                                 << ") ";
  34                         }
  35                         output << endl;
  36                 }
  37                 output << endl;
  38         }
  39         return output;
  40 }
  41
  42
  43 ostream& operator<<(ostream& output, const Interactions::iterator interaction) {
  44         output << (interaction->name.empty() ? "(no name)" : interaction->name)
  45                                         << " [M:" << interaction->multiplicity << ", E:"
  46                                         << interaction->energy << "]\n";
  47         output << endl;
  48         return output;
  49 }
  50
  51
  52 bool Interactions::compareInteraction(const Interaction a, const Interaction b) {
  53         return (a.name < b.name);
  54 }
  55
  56
  57
  58 bool operator==(Directions& firstDirections, Directions& secondDirections) {
  59
  60         /*
  61          * Defining equality operator for Directions:
  62          * A direction is equal to another one, if it contains the same vectors (direction.x/direction.y),
  63          * even if they are arranged in a different order in this direction
  64         */
  65         if (firstDirections.size() != secondDirections.size()) {
  66                 cout << "Directions to compare got different number of vectors" << endl;
  67                 return false;
  68         }
  69
  70         bool equalVector;
  71         int counterEqualVectors = 0;
  72
  73         for (Directions::iterator firstDirection = firstDirections.begin(); firstDirection
  74                         != firstDirections.end(); firstDirection++) {
  75
  76                 for (Directions::iterator secondDirection = secondDirections.begin(); secondDirection
  77                                 != secondDirections.end(); secondDirection++) {
  78                         equalVector = false;
  79                         //cout << "comparing vector 1st direction:(" << firstDirection->x << "/" << firstDirection->y << ") with vector 2nd direction:(" << secondDirection->x << "/" << secondDirection->y << ") " << endl;
  80                         if (firstDirection->x == secondDirection->x && firstDirection->y
  81                                         == secondDirection->y) {
  82                                 equalVector = true;
  83                         }
  84                         if (equalVector) {
  85                                 counterEqualVectors++;
  86                                 break;
  87                         };
  88                 }
  89         }
  90
  91         if (counterEqualVectors == (int) firstDirections.size()) {
  92                 //cout << "Directions equal!" << endl;
  93                 return true;
  94         }
  95         else
  96                 return false;
  97
  98 }
  99
 100
 101 Direction Interaction::mirrorXZ(Directions::iterator direction) {
 102         double newX, newY;
 103         newX = direction->x;
 104         newY = - direction->y;
 105         Direction newDirection(newX, newY);
 106
 107         return newDirection;
 108 }
 109
 110 Direction Interaction::mirrorYZ(Directions::iterator direction) {
 111         double newX, newY;
 112         newX = - direction->x;
 113         newY = direction->y;
 114         Direction newDirection(newX, newY);
 115
 116         return newDirection;
 117 }
 118
 119 vector<Directions> Interaction::rotate90degreeZ(Directions& directions) {
 120
 121         //cout << "Rotations..." << endl;
 122
 123         Directions newDirections;
 124         vector<Directions> returnDirectionsVector;
 125
 126         returnDirectionsVector.push_back(directions);
 127
 128         for (int nrOfRotations = 0; nrOfRotations < 4; nrOfRotations++){
 129
 130                 vector<Directions>::iterator lastRotatedDirections;
 131                 lastRotatedDirections = returnDirectionsVector.end() - 1;
 132
 133                 for (Directions::iterator direction = lastRotatedDirections->begin(); direction != lastRotatedDirections->end(); direction++){
 134                         double newX, newY;
 135                         newX = - direction->y;
 136                         newY = direction->x;
 137                         Direction newDirection(newX, newY);
 138                         newDirections.push_back(newDirection);
 139                 };
 140                 returnDirectionsVector.push_back(newDirections);
 141                 newDirections.clear();
 142         }
 143
 144
 145
 146         return returnDirectionsVector;
 147
 148 }
 149
 150 void Interaction::generateSymmetryEquivalentVecC2v(Directions directions,
 151                 vector<Directions>& symmetryDirections) {
 152
 153         Directions newDirections;
 154         for (Directions::iterator direction = directions.begin(); direction
 155                         != directions.end(); direction++) {
 156                 Direction tempDirection(mirrorXZ(direction));
 157                 //cout << "mirrorXZ yields for: (" << direction->x << "/" << direction->y
 158                 //              << ") the new: (" << tempDirection.x << "/" << tempDirection.y
 159                 //              << ") " << endl;
 160                 newDirections.push_back(tempDirection);
 161         }
 162         symmetryDirections.push_back(newDirections);
 163         newDirections.clear();
 164
 165         for (Directions::iterator direction = directions.begin(); direction
 166                         != directions.end(); direction++) {
 167                 Direction tempDirection(mirrorYZ(direction));
 168                 //cout << "mirrorYZ yields for: (" << direction->x << "/" << direction->y
 169                 //              << ") the new: (" << tempDirection.x << "/" << tempDirection.y
 170                 //              << ") " << endl;
 171                 newDirections.push_back(tempDirection);
 172         }
 173         symmetryDirections.push_back(newDirections);
 174         newDirections.clear();
 175
 176         vector<Directions> newSymmetryDirections;
 177
 178         // Rotating of Symmetry Directions
 179         for (vector<Directions>::iterator symmetryDirection =
 180                         symmetryDirections.begin(); symmetryDirection
 181                         != symmetryDirections.end(); symmetryDirection++) {
 182                         vector<Directions> rotatedDirections;
 183                         rotatedDirections = rotate90degreeZ((*symmetryDirection));
 184                         newSymmetryDirections.insert(newSymmetryDirections.end(), rotatedDirections.begin(), rotatedDirections.end());
 185                         rotatedDirections.clear();
 186         }
 187
 188         // Now rotating the original set of directions finally
 189         vector<Directions> rotatedDirections;
 190         rotatedDirections = rotate90degreeZ(directions);
 191         newSymmetryDirections.insert(newSymmetryDirections.end(), rotatedDirections.begin(), rotatedDirections.end());
 192
 193         // newSymmetryDirections.push_back(rotate90degreeZ(directions));
 194
 195         // Joining newly generated newSymmetryDirections and symmetryDirections
 196         symmetryDirections.insert(symmetryDirections.end(), newSymmetryDirections.begin(), newSymmetryDirections.end());
 197
 198
 199         rotatedDirections.clear();
 200         newSymmetryDirections.clear();
 201
 202         /* DEBUG output
 203         cout << " -------- complete symmetryDirections before checking doubles----- " << endl;
 204         for (vector<Directions>::iterator symmetryVec1 = symmetryDirections.begin(); symmetryVec1
 205                                 != symmetryDirections.end(); symmetryVec1++) {
 206                 //cout << "new interaction" << " with " << symmetryVec1->size() << " as size" << endl;
 207                 for (Directions::iterator direction = symmetryVec1->begin(); direction != symmetryVec1->end(); direction++){
 208                         //cout << ". - vector (" << direction->x << "/" << direction->y << ")" << endl;
 209                 }
 210         }
 211
 212         cout << " -------- complete symmetryDirections  end ------------------------ " << endl << endl;
 213         */
 214         // sorting out double occurrences of interaction vectors
 215         //cout << "sorting out double occurencies!" << endl << endl << endl;
 216
 217
 218         //cout << symmetryDirections.size() << " symmetrized Directions after rotations and mirroring" << endl;
 219
 220
 221         vector<Directions> uniqueSymmetryDirections;
 222
 223         for (vector<Directions>::iterator symmetryDirections1 =
 224                         symmetryDirections.begin(); symmetryDirections1
 225                         != symmetryDirections.end(); symmetryDirections1++) {
 226                 //cout << "comparing to original directions" << endl;
 227                 if (*symmetryDirections1 == directions) {
 228                         continue;
 229                         };
 230                 if (uniqueSymmetryDirections.empty())
 231                         uniqueSymmetryDirections.push_back(*symmetryDirections1);
 232
 233                 vector<Directions> comparisonDirections;
 234                 comparisonDirections.insert(comparisonDirections.end(), uniqueSymmetryDirections.begin(), uniqueSymmetryDirections.end());
 235                 bool alreadyPresentInteraction = false;
 236                 for (vector<Directions>::iterator symmetryDirections2 =
 237                         comparisonDirections.begin(); symmetryDirections2
 238                         != comparisonDirections.end(); symmetryDirections2++) {
 239                                 //cout << " comparing symmetryDirections " << endl;
 240                                 //cout << " symmetryDirections2 size: " << symmetryDirections2->size() << " with symmetryDirections1 size: " << symmetryDirections1->size() << endl;
 241                                 if ( (*symmetryDirections2 == *symmetryDirections1) ) {
 242                                         alreadyPresentInteraction = true;
 243                                 };
 244                 };
 245                 if (!alreadyPresentInteraction){
 246                                 //cout << " pushing back direction " << endl;
 247                                 uniqueSymmetryDirections.push_back(*symmetryDirections1);
 248                 };
 249
 250                 comparisonDirections.clear();
 251
 252         };
 253         //cout << "finished sorting out" << endl;
 254
 255         // assigning sorted directions!
 256         symmetryDirections = uniqueSymmetryDirections;
 257
 258         /* DEBUG
 259         cout << " -------- complete symmetryDirections after checking doubles----- " << endl;
 260
 261         for (vector<Directions>::iterator symmetryVec1 = symmetryDirections.begin(); symmetryVec1
 262                                 != symmetryDirections.end(); symmetryVec1++) {
 263                 cout << "new interaction" << endl;
 264                 for (Directions::iterator direction = symmetryVec1->begin(); direction != symmetryVec1->end(); direction++){
 265                         cout << " - vector (" << direction->x << "/" << direction->y << ")" << endl;
 266                 }
 267         }
 268
 269         cout << " -------- complete symmetryDirections after doubles checking end ------------------------ " << endl << endl;
 270         */
 271
 272         return;
 273 }
 274
 275
 276 void Interactions::symmetrizeInteractions(Interactions& interactions) {
 277         //cout << "#############     calling symmetrizer  #################" << endl;
 278         for (Interactions::iterator interaction = interactions.begin(); interaction != interactions.end(); interaction++){
 279                 //cout << "working on interaction " << interaction->name << endl;
 280                 interaction->generateSymmetryEquivalentVecC2v(interaction->directions, interaction->symmetryDirections);
 281         };
 282
 283         return;
 284 }
 285
 286 void Interactions::joinSymmetryEquivalent(Interactions& interactions) {
 287         /* joining the equivalent representations of multi body interactions into one interaction object
 288          * we are doing this, by going through interactions comparing the name
 289          * if it is the same name, we are joining the direction and symmetryDirections object of the second one into the
 290          * symmetryDirections object of the first one
 291          * then the second direction gets deleted
 292          *
 293          */
 294         //cout << " -------- joining equivalent representations! --------" << endl;
 295
 296
 297         /*
 298         Interactions::iterator interaction = interactions.begin();
 299         Interactions joinedInteractions;
 300         while (interaction != interactions.end()) {
 301                 if (joinedInteractions.empty()) {
 302                         joinedInteractions.push_back(*interaction);
 303                         interaction++;
 304                         continue;
 305                 } else { // checking for members of joinedInteractions...
 306                         Interactions comparisonInteractions;
 307                         comparisonInteractions.clear();
 308                         comparisonInteractions.insert(comparisonInteractions.begin(),
 309                                         joinedInteractions.begin(), joinedInteractions.end());
 310                         for (Interactions::iterator comparisonInteraction =
 311                                         comparisonInteractions.begin(); comparisonInteraction
 312                                         != comparisonInteractions.end(); comparisonInteraction++) {
 313                                 cout << "interaction is: " << interaction->name << endl;
 314                                 cout << "comparing to interaction name: " << comparisonInteraction->name << endl;
 315                                 if (comparisonInteraction->name.compare(interaction->name) == 0 // exactly equal
 316                                                 && comparisonInteraction->directions // nothing to do
 317                                                                 == interaction->directions)
 318                                         continue;
 319                                 else if (comparisonInteraction->name.compare(interaction->name)
 320                                                 == 0) { // equal representations - joining necessary
 321                                         Interactions::iterator joinedInteraction =
 322                                                         joinedInteractions.begin();
 323                                         while (joinedInteraction->name.compare(
 324                                                         comparisonInteraction->name) != 0)
 325                                                 joinedInteraction++;
 326                                         cout << "joining equivalent representations" << endl;
 327                                         joinedInteraction->symmetryDirections.push_back(
 328                                                         interaction->directions);
 329                                         cout << "pushed back original directions" << endl;
 330                                         joinedInteraction->symmetryDirections.insert(
 331                                                         joinedInteraction->symmetryDirections.end(),
 332                                                         interaction->symmetryDirections.begin(),
 333                                                         interaction->symmetryDirections.end());
 334                                         cout << "pushed back symmetryDirections directions" << endl;
 335                                         //if ( interaction != interactions.end())
 336                                         //                      interaction++;
 337                                         continue;
 338                                         }
 339                                 };
 340                         }
 341                 if (interaction != interactions.end()){
 342                         joinedInteractions.push_back(*interaction);
 343                         cout << "pushed back original interaction to joinedInteractions" << endl;
 344                         interaction++;
 345                 };
 346         };
 347         */
 348
 349         // new version of joining!
 350
 351         Interactions joinedInteractions;
 352         if (joinedInteractions.empty())
 353                         joinedInteractions.push_back(*(interactions.begin()));
 354
 355         Interactions::iterator joinedInteraction = joinedInteractions.begin();
 356         Interactions::iterator joinedInteraction_end = joinedInteractions.end();
 357         int maxChecks = interactions.size();
 358         int nrOfChecks = 0;
 359
 360         map<string, bool> pushedInteractions;
 361
 362         pushedInteractions[interactions.begin()->name] = true;
 363
 364
 365
 366
 367         while (joinedInteraction != joinedInteraction_end && nrOfChecks <= maxChecks) {
 368                 string workingInteractionName = joinedInteraction->name;
 369                 for (Interactions::iterator interaction = interactions.begin(); interaction != interactions.end(); interaction++) {
 370                         //cout << "joinedInteraction is: " << joinedInteraction->name << "\n";
 371                         //cout << "Interaction is      : " << interaction->name << "\n";
 372                         if (joinedInteraction->name.compare(interaction->name) == 0 && joinedInteraction->directions == interaction->directions)
 373                         continue;
 374                         if (joinedInteraction->name.compare(interaction->name) == 0) {
 375                                 joinedInteraction->symmetryDirections.push_back(
 376                                                 interaction->directions);
 377                                 //cout << "pushed back original directions" << endl;
 378                                 joinedInteraction->symmetryDirections.insert(
 379                                                 joinedInteraction->symmetryDirections.end(),
 380                                                 interaction->symmetryDirections.begin(),
 381                                                 interaction->symmetryDirections.end());
 382                                 //cout << "pushed back symmetryDirections directions" << endl;
 383                                 continue;
 384                         }
 385                         else { // checking if already pushed to joinedInteractions
 386                                 if ((pushedInteractions[interaction->name]))
 387                                         continue;
 388                                 else {
 389                                         //cout << " pushing back interaction : " << interaction->name << " to joinedInteractions \n";
 390                                         pushedInteractions[interaction->name] = true;
 391                                         joinedInteractions.push_back(*interaction); // push_back destroys all iterators
 392                                         joinedInteraction = joinedInteractions.begin(); // retrieveing valid iterators
 393                                         joinedInteraction_end = joinedInteractions.end();
 394                                         while (joinedInteraction->name.compare(workingInteractionName) != 0)
 395                                                 joinedInteraction++; // now being at the original position again
 396                                         // joinedInteraction++; // new position for comparing
 397                                         continue;
 398                                 };
 399                         };
 400                 };
 401                 joinedInteraction++;
 402                 nrOfChecks++;
 403         };
 404
 405
 406         // interactions.assign(joinedInteractions.begin(), joinedInteractions.end());
 407         //cout << "interactions before consolidating" << endl;
 408         //cout << interactions << endl;
 409
 410         //cout << "interactions after consolidating" << endl;
 411         //cout << joinedInteractions << endl;
 412         interactions = joinedInteractions;
 413
 414         // new!
 415         //cout << "now checking once more for equivalent directions in symmetryDirections! \n";
 416
 417         for (Interactions::iterator interaction = interactions.begin(); interaction
 418                         != interactions.end(); interaction++) {
 419                 //cout << "working on interaction " << interaction->name << endl;
 420                 Directions directions = interaction->directions;
 421                 vector<Directions> symmetryDirections = interaction->symmetryDirections;
 422
 423                 vector<Directions> uniqueSymmetryDirections;
 424
 425                 for (vector<Directions>::iterator symmetryDirections1 =
 426                                 symmetryDirections.begin(); symmetryDirections1
 427                                 != symmetryDirections.end(); symmetryDirections1++) {
 428                         //cout << "comparing to original directions" << endl;
 429                         if (*symmetryDirections1 == directions) {
 430                                 continue;
 431                         };
 432                         if (uniqueSymmetryDirections.empty())
 433                                 uniqueSymmetryDirections.push_back(*symmetryDirections1);
 434
 435                         vector<Directions> comparisonDirections;
 436                         comparisonDirections.insert(comparisonDirections.end(),
 437                                         uniqueSymmetryDirections.begin(),
 438                                         uniqueSymmetryDirections.end());
 439                         bool alreadyPresentInteraction = false;
 440                         for (vector<Directions>::iterator symmetryDirections2 =
 441                                         comparisonDirections.begin(); symmetryDirections2
 442                                         != comparisonDirections.end(); symmetryDirections2++) {
 443                                 //cout << " comparing symmetryDirections " << endl;
 444                                 //cout << " symmetryDirections2 size: "
 445                                 //              << symmetryDirections2->size()
 446                                 //              << " with symmetryDirections1 size: "
 447                                 //              << symmetryDirections1->size() << endl;
 448                                 if ((*symmetryDirections2 == *symmetryDirections1)) {
 449                                         alreadyPresentInteraction = true;
 450                                 };
 451                         };
 452                         if (!alreadyPresentInteraction) {
 453                                 //cout << " pushing back direction " << endl;
 454                                 uniqueSymmetryDirections.push_back(*symmetryDirections1);
 455                         };
 456
 457                         comparisonDirections.clear();
 458
 459                 };
 460                 // assigning sorted directions!
 461                 interaction->symmetryDirections = uniqueSymmetryDirections;
 462         };
 463
 464
 465
 466 }
 467
 468 void Interactions::generateCombinations(){
 469         /* this routine will use next_combination() from combination.hpp to generate all possible combinations of
 470          * a given size out of the defined interations
 471          * the purpose is to find the "best" set of interactions that is giving the lowest CV Score
 472          * to describe the total energies of the input structures
 473          *
 474          * the original routine is from Ben Bear, Herve Bronnimann (2007)
 475          * http://photon.poly.edu/~hbr/boost/combination.hpp
 476          */
 477
 478         const int r = 48; // size for the combination to be build = will start from 1 and will go to interactions.size()-1
 479         const int n = 48; // total number of elements in basis set = total number of interactions
 480          std::vector<int> v_int(n);
 481          for (int i = 0; i < n; ++i) { v_int[i] = i; } // generating starting vector for combinations
 482
 483          int N = 0;
 484          do {
 485              ++N;
 486              std::cout << "[ " << v_int[0];
 487                  for (int j = 1; j < r; ++j) { std::cout << ", " << v_int[j]; }
 488                  std::cout << " ]" << std::endl;
 489
 490          } while (next_combination(v_int.begin(), v_int.begin() + r, v_int.end()));
 491          std::cout << "Found " << N << " combinations of size " << r << " without repetitions"
 492                    << " from a set of " << n << " elements." << std::endl;
 493 }
 494