Line data    Source code

       1             : /*
       2             :  * Copyright (c) 2016 Sebastian Weber, Henri Menke. All rights reserved.
       3             :  *
       4             :  * This file is part of the pairinteraction library.
       5             :  *
       6             :  * The pairinteraction library is free software: you can redistribute it and/or modify
       7             :  * it under the terms of the GNU Lesser General Public License as published by
       8             :  * the Free Software Foundation, either version 3 of the License, or
       9             :  * (at your option) any later version.
      10             :  *
      11             :  * The pairinteraction library is distributed in the hope that it will be useful,
      12             :  * but WITHOUT ANY WARRANTY; without even the implied warranty of
      13             :  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
      14             :  * GNU Lesser General Public License for more details.
      15             :  *
      16             :  * You should have received a copy of the GNU Lesser General Public License
      17             :  * along with the pairinteraction library. If not, see <http://www.gnu.org/licenses/>.
      18             :  */
      19             : 
      20             : #include "Basisnames.hpp"
      21             : #include <stdexcept>
      22             : 
      23           7 : BasisnamesOne::BasisnamesOne() = default;
      24           0 : BasisnamesOne BasisnamesOne::fromStates(const std::vector<StateOneOld> &names) {
      25           0 :     BasisnamesOne basisnames;
      26           0 :     basisnames.names_ = names;
      27           0 :     basisnames.dim_ = names.size();
      28           0 :     return basisnames;
      29             : }
      30           0 : BasisnamesOne BasisnamesOne::fromFirst(const Configuration &config) {
      31           0 :     StateOneOld startstate;
      32           0 :     config["n1"] >> startstate.n;
      33           0 :     config["l1"] >> startstate.l;
      34           0 :     config["j1"] >> startstate.j;
      35           0 :     config["m1"] >> startstate.m;
      36             : 
      37           0 :     BasisnamesOne basisnames;
      38           0 :     basisnames._constructedFromFirst = true;
      39           0 :     basisnames.configure(config);
      40           0 :     basisnames.build(startstate, config["species1"].str());
      41           0 :     return basisnames;
      42             : }
      43           2 : BasisnamesOne BasisnamesOne::fromFirst(const std::shared_ptr<const BasisnamesTwo> &basis_two) {
      44           4 :     Configuration config = basis_two->getConf();
      45           4 :     StateOneOld startstate;
      46           2 :     config["n1"] >> startstate.n;
      47           2 :     config["l1"] >> startstate.l;
      48           2 :     config["j1"] >> startstate.j;
      49           2 :     config["m1"] >> startstate.m;
      50             : 
      51           2 :     BasisnamesOne basisnames;
      52           2 :     basisnames._constructedFromFirst = true;
      53           2 :     basisnames.configure(config);
      54           2 :     basisnames.build(startstate, config["species1"].str(), basis_two, 0);
      55           4 :     return basisnames;
      56             : }
      57           0 : BasisnamesOne BasisnamesOne::fromSecond(const Configuration &config) {
      58           0 :     StateOneOld startstate;
      59           0 :     config["n2"] >> startstate.n;
      60           0 :     config["l2"] >> startstate.l;
      61           0 :     config["j2"] >> startstate.j;
      62           0 :     config["m2"] >> startstate.m;
      63             : 
      64           0 :     BasisnamesOne basisnames;
      65           0 :     basisnames._constructedFromFirst = false;
      66           0 :     basisnames.configure(config);
      67           0 :     basisnames.build(startstate, config["species2"].str());
      68           0 :     return basisnames;
      69             : }
      70           2 : BasisnamesOne BasisnamesOne::fromSecond(const std::shared_ptr<const BasisnamesTwo> &basis_two) {
      71           4 :     Configuration config = basis_two->getConf();
      72           4 :     StateOneOld startstate;
      73           2 :     config["n2"] >> startstate.n;
      74           2 :     config["l2"] >> startstate.l;
      75           2 :     config["j2"] >> startstate.j;
      76           2 :     config["m2"] >> startstate.m;
      77             : 
      78           2 :     BasisnamesOne basisnames;
      79           2 :     basisnames._constructedFromFirst = false;
      80           2 :     basisnames.configure(config);
      81           2 :     basisnames.build(startstate, config["species2"].str(), basis_two, 1);
      82           4 :     return basisnames;
      83             : }
      84           3 : BasisnamesOne BasisnamesOne::fromBoth(const Configuration &config) {
      85           6 :     StateTwoOld startstate;
      86           3 :     config["n1"] >> startstate.n[0];
      87           3 :     config["l1"] >> startstate.l[0];
      88           3 :     config["j1"] >> startstate.j[0];
      89           3 :     config["m1"] >> startstate.m[0];
      90           3 :     config["n2"] >> startstate.n[1];
      91           3 :     config["l2"] >> startstate.l[1];
      92           3 :     config["j2"] >> startstate.j[1];
      93           3 :     config["m2"] >> startstate.m[1];
      94             : 
      95           3 :     if (config["species1"].str() != config["species2"].str()) {
      96           0 :         throw std::runtime_error(
      97           0 :             "BasisnamesOne::fromBoth can only be used if both atoms are of the same species.");
      98             :     }
      99             : 
     100           3 :     BasisnamesOne basisnames;
     101           3 :     basisnames._constructedFromFirst = false;
     102           3 :     basisnames.configure(config);
     103           3 :     basisnames.build(startstate.order(), config["species1"].str());
     104           6 :     return basisnames;
     105             : 
     106             :     /*if ((startstate.n[0] == startstate.n[1]) &&
     107             :         (startstate.l[0] == startstate.l[1]) &&
     108             :         (startstate.j[0] == startstate.j[1]) &&
     109             :         (startstate.m[0] == startstate.m[1])) {
     110             :         BasisnamesOne basisnames;
     111             :         basisnames._constructedFromFirst = false;
     112             :         basisnames.configure(config);
     113             :         basisnames.build(startstate.first(), config["species1"].str());
     114             :         return basisnames;
     115             :     } else {
     116             : 
     117             :     }*/ // TODO
     118             : }
     119           3 : void BasisnamesOne::build(StateTwoOld startstate, const std::string &species) {
     120           3 :     states_initial.push_back(startstate.first());  // TODO correct for idx
     121           3 :     states_initial.push_back(startstate.second()); // TODO correct for idx
     122             : 
     123           3 :     conf["species1"] << species;
     124           3 :     conf["n1"] << startstate.n[0];
     125           3 :     conf["l1"] << startstate.l[0];
     126           3 :     conf["j1"] << startstate.j[0];
     127           3 :     conf["m1"] << startstate.m[0];
     128           3 :     conf["n2"] << startstate.n[1];
     129           3 :     conf["l2"] << startstate.l[1];
     130           3 :     conf["j2"] << startstate.j[1];
     131           3 :     conf["m2"] << startstate.m[1];
     132             : 
     133           6 :     std::unordered_set<StateOneOld> names_set;
     134             : 
     135           3 :     idx_t idx = 0;
     136             : 
     137           3 :     if (delta_l < 0) {
     138           0 :         delta_l = std::fmax(startstate.l[0], startstate.l[1]) +
     139           0 :             std::fmax(startstate.n[0], startstate.n[1]) + delta_n - 1;
     140             :     }
     141           3 :     if (delta_j < 0) {
     142           3 :         delta_j = std::fmax(startstate.j[0], startstate.j[1]) +
     143           3 :             std::fmax(startstate.n[0], startstate.n[1]) + delta_n - 0.5;
     144             :     }
     145           3 :     if (delta_m < 0) {
     146           3 :         delta_m = std::fmax(startstate.m[0], startstate.m[1]) +
     147           3 :             std::fmax(startstate.n[0], startstate.n[1]) + delta_n - 0.5;
     148             :     }
     149             : 
     150             :     // loop over quantum numbers of startstate1
     151          24 :     for (int n = std::fmax(0, startstate.n[0] - delta_n); n <= startstate.n[0] + delta_n; ++n) {
     152         126 :         for (int l = std::fmax(0, startstate.l[0] - delta_l);
     153         126 :              l <= std::fmin(n - 1, startstate.l[0] + delta_l); ++l) {
     154         294 :             for (float j = std::fmax(std::fabs(l - 0.5), startstate.j[0] - delta_j);
     155         294 :                  j <= std::fmin(l + 0.5, startstate.j[0] + delta_j); ++j) {
     156        1239 :                 for (float m = std::fmax(-j, startstate.m[0] - delta_m);
     157        1239 :                      m <= std::fmin(j, startstate.m[0] + delta_m); ++m) {
     158        1050 :                     auto result = names_set.insert(StateOneOld(idx, n, l, j, m));
     159        1050 :                     if (result.second) {
     160        1050 :                         idx++;
     161             :                     }
     162             :                 }
     163             :             }
     164             :         }
     165             :     }
     166             : 
     167             :     // loop over quantum numbers of startstate2
     168          24 :     for (int n = std::fmax(0, startstate.n[1] - delta_n); n <= startstate.n[1] + delta_n; ++n) {
     169         105 :         for (int l = std::fmax(0, startstate.l[1] - delta_l);
     170         105 :              l <= std::fmin(n - 1, startstate.l[1] + delta_l); ++l) {
     171         231 :             for (float j = std::fmax(std::fabs(l - 0.5), startstate.j[1] - delta_j);
     172         231 :                  j <= std::fmin(l + 0.5, startstate.j[1] + delta_j); ++j) {
     173         819 :                 for (float m = std::fmax(-j, startstate.m[1] - delta_m);
     174         819 :                      m <= std::fmin(j, startstate.m[1] + delta_m); ++m) {
     175         672 :                     auto result = names_set.insert(StateOneOld(idx, n, l, j, m));
     176         672 :                     if (result.second) {
     177          96 :                         idx++;
     178             :                     }
     179             :                 }
     180             :             }
     181             :         }
     182             :     }
     183             : 
     184           3 :     std::set<StateOneOld> names_ordered(names_set.begin(), names_set.end());
     185           3 :     names_ = std::vector<StateOneOld>(names_ordered.begin(), names_ordered.end());
     186             : 
     187           3 :     dim_ = idx;
     188           3 : }
     189           0 : void BasisnamesOne::build(StateOneOld startstate, const std::string &species) {
     190           0 :     states_initial.push_back(startstate); // TODO correct for idx
     191             : 
     192           0 :     conf["species1"] << species;
     193           0 :     conf["n1"] << startstate.n;
     194           0 :     conf["l1"] << startstate.l;
     195           0 :     conf["j1"] << startstate.j;
     196           0 :     conf["m1"] << startstate.m;
     197           0 :     conf["n2"] << "";
     198           0 :     conf["l2"] << "";
     199           0 :     conf["j2"] << "";
     200           0 :     conf["m2"] << "";
     201             : 
     202           0 :     idx_t idx = 0;
     203             : 
     204           0 :     if (delta_l < 0) {
     205           0 :         delta_l = startstate.l + startstate.n + delta_n - 1;
     206             :     }
     207           0 :     if (delta_j < 0) {
     208           0 :         delta_j = startstate.j + startstate.n + delta_n - 0.5;
     209             :     }
     210           0 :     if (delta_m < 0) {
     211           0 :         delta_m = startstate.m + startstate.n + delta_n - 0.5;
     212             :     }
     213             : 
     214             :     // loop over quantum numbers
     215           0 :     for (int n = std::fmax(0, startstate.n - delta_n); n <= startstate.n + delta_n; ++n) {
     216           0 :         for (int l = std::fmax(0, startstate.l - delta_l);
     217           0 :              l <= std::fmin(n - 1, startstate.l + delta_l); ++l) {
     218           0 :             for (float j = std::fmax(std::fabs(l - 0.5), startstate.j - delta_j);
     219           0 :                  j <= std::fmin(l + 0.5, startstate.j + delta_j); ++j) {
     220           0 :                 for (float m = std::fmax(-j, startstate.m - delta_m);
     221           0 :                      m <= std::fmin(j, startstate.m + delta_m); ++m) { // TODO
     222           0 :                     names_.emplace_back(idx++, n, l, j, m);
     223             :                 }
     224             :             }
     225             :         }
     226             :     }
     227             : 
     228           0 :     dim_ = idx;
     229           0 : }
     230           4 : void BasisnamesOne::build(StateOneOld startstate, const std::string &species,
     231             :                           const std::shared_ptr<const BasisnamesTwo> &basis_two, int i) {
     232           4 :     states_initial.push_back(startstate); // TODO correct for idx
     233             : 
     234           4 :     conf["species1"] << species;
     235           4 :     conf["n1"] << startstate.n;
     236           4 :     conf["l1"] << startstate.l;
     237           4 :     conf["j1"] << startstate.j;
     238           4 :     conf["m1"] << startstate.m;
     239           4 :     conf["n2"] << "";
     240           4 :     conf["l2"] << "";
     241           4 :     conf["j2"] << "";
     242           4 :     conf["m2"] << "";
     243             : 
     244           8 :     std::unordered_set<StateOneOld> names_set;
     245             : 
     246           4 :     idx_t idx = 0;
     247             : 
     248             :     // loop over quantum numbers
     249      100356 :     for (auto state : *basis_two) {
     250             :         auto result =
     251      100352 :             names_set.insert(StateOneOld(idx, state.n[i], state.l[i], state.j[i], state.m[i]));
     252      100352 :         if (result.second) {
     253         448 :             idx++;
     254             :         }
     255             :     }
     256             : 
     257           4 :     std::set<StateOneOld> names_ordered(names_set.begin(), names_set.end());
     258           4 :     names_ = std::vector<StateOneOld>(names_ordered.begin(), names_ordered.end());
     259             : 
     260           4 :     dim_ = idx;
     261           4 : }
     262           8 : const std::vector<StateOneOld> &BasisnamesOne::initial() const { return states_initial; }
     263           5 : void BasisnamesOne::removeUnnecessaryStates(const std::vector<bool> &is_necessary) {
     264          10 :     auto tmp = names_;
     265           5 :     names_.clear();
     266           5 :     names_.reserve(tmp.size());
     267             : 
     268             :     // loop over all one-atom states
     269           5 :     idx_t idx = 0;
     270        1375 :     for (auto state : tmp) {
     271        1370 :         if (is_necessary[state.idx]) {
     272         598 :             state.idx = idx;
     273         598 :             names_.push_back(state);
     274             :             // TODO update indices of states_initial
     275         598 :             ++idx;
     276             :         }
     277             :     }
     278             : 
     279           5 :     dim_ = idx;
     280           5 :     names_.shrink_to_fit();
     281           5 : }
     282             : 
     283           0 : bool BasisnamesOne::constructedFromFirst() { return _constructedFromFirst; }
     284             : 
     285           3 : void BasisnamesOne::save(const std::string &path) {
     286           6 :     std::ofstream csvfile;
     287           3 :     csvfile.open(path);
     288         377 :     for (const auto &state : *this) {
     289         374 :         csvfile << state.idx << "\t" << state.n << "\t" << state.l << "\t" << state.j << "\t"
     290         374 :                 << state.m << std::endl;
     291             :     }
     292           3 :     csvfile.close();
     293           3 : }
     294             : 
     295           4 : BasisnamesTwo::BasisnamesTwo(const std::shared_ptr<const BasisnamesOne> &basis_one1) {
     296           8 :     const Configuration conf1 = basis_one1->getConf();
     297             : 
     298           4 :     if (conf1["n2"].str().empty()) {
     299           0 :         throw std::runtime_error("BasisnamesTwo can be only constructed from two "
     300           0 :                                  "BasisnamesOne::fromFirst / BasisnamesOne::fromSecond.");
     301             :     }
     302             : 
     303           4 :     configure(conf1);
     304           4 :     conf["combined"] << 1;
     305             : 
     306           8 :     StateTwoOld startstate;
     307           4 :     conf1["n1"] >> startstate.n[0];
     308           4 :     conf1["l1"] >> startstate.l[0];
     309           4 :     conf1["j1"] >> startstate.j[0];
     310           4 :     conf1["m1"] >> startstate.m[0];
     311           4 :     conf1["n2"] >> startstate.n[1];
     312           4 :     conf1["l2"] >> startstate.l[1];
     313           4 :     conf1["j2"] >> startstate.j[1];
     314           4 :     conf1["m2"] >> startstate.m[1];
     315             : 
     316             :     std::array<std::string, 2> species(
     317           8 :         {{conf1["species1"].str(),
     318          12 :           conf1["species1"].str()}}); // TODO : species in state class mit aufnehmen
     319           4 :     build(startstate, species, basis_one1, basis_one1);
     320           4 : }
     321             : 
     322           0 : BasisnamesTwo::BasisnamesTwo(const std::shared_ptr<const BasisnamesOne> &basis_one1,
     323           0 :                              const std::shared_ptr<const BasisnamesOne> &basis_one2) {
     324           0 :     const Configuration conf1 = basis_one1->getConf();
     325           0 :     const Configuration conf2 = basis_one2->getConf();
     326             : 
     327           0 :     if (!conf1["n2"].str().empty() || !conf2["n2"].str().empty()) {
     328           0 :         throw std::runtime_error(
     329           0 :             "BasisnamesTwo can be only constructed from one single BasisnamesOne::fromBoth.");
     330             :     }
     331             : 
     332           0 :     configure(conf1);
     333           0 :     conf["combined"] << 0;
     334             : 
     335           0 :     StateTwoOld startstate;
     336           0 :     conf1["n1"] >> startstate.n[0];
     337           0 :     conf1["l1"] >> startstate.l[0];
     338           0 :     conf1["j1"] >> startstate.j[0];
     339           0 :     conf1["m1"] >> startstate.m[0];
     340           0 :     conf2["n1"] >> startstate.n[1];
     341           0 :     conf2["l1"] >> startstate.l[1];
     342           0 :     conf2["j1"] >> startstate.j[1];
     343           0 :     conf2["m1"] >> startstate.m[1];
     344             : 
     345             :     std::array<std::string, 2> species(
     346           0 :         {{conf1["species1"].str(),
     347           0 :           conf2["species1"].str()}}); // TODO : species in state class mit aufnehmen
     348           0 :     build(startstate, species, basis_one1, basis_one2);
     349           0 : }
     350             : 
     351           2 : const StateTwoOld &BasisnamesTwo::initial() const { return state_initial; }
     352             : 
     353           0 : void BasisnamesTwo::removeUnnecessaryStates(const std::vector<bool> &is_necessary) {
     354           0 :     auto tmp = names_;
     355           0 :     names_.clear();
     356           0 :     names_.reserve(tmp.size());
     357             : 
     358             :     // loop over all two-atom states
     359           0 :     bool state_initial_found = false;
     360           0 :     idx_t idx = 0;
     361           0 :     for (auto state : tmp) {
     362           0 :         if (is_necessary[state.idx]) {
     363           0 :             state.idx = idx;
     364           0 :             names_.push_back(state);
     365             : 
     366           0 :             if (!state_initial_found && state == state_initial) {
     367           0 :                 state_initial.idx = idx;
     368           0 :                 state_initial_found = true;
     369             :             }
     370             : 
     371           0 :             ++idx;
     372             :         }
     373             :     }
     374             : 
     375           0 :     dim_ = idx;
     376           0 :     names_.shrink_to_fit();
     377           0 : }
     378             : 
     379           0 : void BasisnamesTwo::removeUnnecessaryStatesKeepIdx(const std::vector<bool> &is_necessary) {
     380           0 :     auto tmp = names_;
     381           0 :     names_.clear();
     382           0 :     names_.reserve(tmp.size());
     383             : 
     384             :     // loop over all two-atom states
     385           0 :     for (auto state : tmp) {
     386           0 :         if (is_necessary[state.idx]) {
     387           0 :             names_.push_back(state);
     388             :         }
     389             :     }
     390           0 :     names_.shrink_to_fit();
     391           0 : }
     392             : 
     393           4 : void BasisnamesTwo::build(StateTwoOld startstate, std::array<std::string, 2> species,
     394             :                           const std::shared_ptr<const BasisnamesOne> &basis_one1,
     395             :                           const std::shared_ptr<const BasisnamesOne> &basis_one2) {
     396           4 :     state_initial = startstate;
     397             : 
     398           4 :     conf["species1"] << species[0];
     399           4 :     conf["n1"] << startstate.n[0];
     400           4 :     conf["l1"] << startstate.l[0];
     401           4 :     conf["j1"] << startstate.j[0];
     402           4 :     conf["m1"] << startstate.m[0];
     403           4 :     conf["species2"] << species[1];
     404           4 :     conf["n2"] << startstate.n[1];
     405           4 :     conf["l2"] << startstate.l[1];
     406           4 :     conf["j2"] << startstate.j[1];
     407           4 :     conf["m2"] << startstate.m[1];
     408             : 
     409           4 :     size_t size = basis_one1->size() * basis_one2->size();
     410           4 :     names_.reserve(size);
     411             : 
     412           4 :     idx_t idx = 0;
     413             : 
     414             :     // loop over single atom states
     415           4 :     bool state_initial_found = false;
     416         452 :     for (const auto &state_1 : *basis_one1) {
     417      100800 :         for (const auto &state_2 : *basis_one2) {
     418      100352 :             names_.emplace_back(idx, state_1, state_2);
     419             : 
     420      100352 :             if (!state_initial_found && names_.back() == state_initial) {
     421           2 :                 state_initial.idx = idx;
     422           2 :                 state_initial_found = true;
     423             :             }
     424             : 
     425      100352 :             ++idx;
     426             :         }
     427             :     }
     428             : 
     429           4 :     dim_ = idx;
     430           4 : }
     431             : 
     432           2 : void BasisnamesTwo::save(const std::string &path) {
     433           4 :     std::ofstream csvfile;
     434           2 :     csvfile.open(path);
     435       50178 :     for (const auto &state : *this) {
     436       50176 :         csvfile << state.idx << "\t" << state.n[0] << "\t" << state.l[0] << "\t" << state.j[0]
     437       50176 :                 << "\t" << state.m[0] << "\t" << state.n[1] << "\t" << state.l[1] << "\t"
     438       50176 :                 << state.j[1] << "\t" << state.m[1] << std::endl;
     439             :     }
     440           2 :     csvfile.close();
     441           2 : }