1 // 2 // BAGEL - Brilliantly Advanced General Electronic Structure Library 3 // Filename: quatmatrix.h 4 // Copyright (C) 2014 Toru Shiozaki 5 // 6 // Author: Toru Shiozaki <shiozaki@northwestern.edu> 7 // Maintainer: Shiozaki group 8 // 9 // This file is part of the BAGEL package. 10 // 11 // This program is free software: you can redistribute it and/or modify 12 // it under the terms of the GNU General Public License as published by 13 // the Free Software Foundation, either version 3 of the License, or 14 // (at your option) any later version. 15 // 16 // This program is distributed in the hope that it will be useful, 17 // but WITHOUT ANY WARRANTY; without even the implied warranty of 18 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 19 // GNU General Public License for more details. 20 // 21 // You should have received a copy of the GNU General Public License 22 // along with this program. If not, see <http://www.gnu.org/licenses/>. 23 // 24 // 25 #ifndef __SRC_MATH_QUATMATRIX_H 26 #define __SRC_MATH_QUATMATRIX_H 27 28 #include <src/util/math/zmatrix.h> 29 #include <src/util/math/algo.h> 30 31 namespace bagel { 32 33 class QuatMatrix : public ZMatrix { 34 public: QuatMatrix(const ZMatrix & o)35 QuatMatrix(const ZMatrix& o) : ZMatrix(o) { } QuatMatrix(ZMatrix && o)36 QuatMatrix(ZMatrix&& o) : ZMatrix(std::move(o)) { } 37 // TODO : implement constructor that can take "00" and "01" matrices and build the rest of the matrix 38 diagonalize(VecView eig)39 void diagonalize(VecView eig) override { 40 #ifndef NDEBUG 41 auto tmp = std::make_shared<ZMatrix>(*this); 42 VectorB eig2(ndim()); 43 std::static_pointer_cast<ZMatrix>(tmp)->diagonalize(eig2); 44 #endif 45 assert(ndim() == mdim()); 46 assert(eig.size() >= ndim()); 47 // assert that matrix is hermitian to ensure real eigenvalues 48 assert(is_hermitian(1.0e-10)); 49 50 zquatev(ndim(), data(), ndim(), eig.data()); 51 52 // zquatev_ only gives half the eigenvalues; get the others using symmetry 53 for (int i = 0; i != ndim()/2; ++i) { 54 eig(ndim()/2+i) = eig(i); 55 #ifndef NDEBUG 56 if (std::max(std::abs(eig(i)-eig2(i*2)), std::abs(eig(i)-eig2(i*2+1))) > 1.0e-6) 57 std::cout << " warning - eigenvalues between quaternion and standard diagonalization do not match:" << std::setprecision(10) << std::fixed << std::setw(20) << eig(i) << std::setw(20) << eig2(i*2) << std::setw(20) << eig2(i*2+1) << std::endl; 58 #endif 59 } 60 61 synchronize(); 62 } 63 synchronize()64 void synchronize() { 65 mpi__->broadcast(data(), size(), 0); 66 } 67 68 // Average out errors in time-reversal symmetry t_symmetrize()69 void t_symmetrize() { 70 assert(mdim()%2 == 0 && ndim()%2 == 0); 71 const size_t m = mdim()/2; 72 const size_t n = ndim()/2; 73 for (size_t i=0; i!=n; ++i) { 74 for (size_t j=0; j!=m; ++j) { 75 element(i, j) = 0.5*(element(i, j)+std::conj(element(n+i, m+j))); 76 element(n+i, j) = 0.5*(element(n+i, j)-std::conj(element(i, m+j))); 77 element(n+i, m+j) = std::conj(element(i, j)); 78 element(i, m+j) = -std::conj(element(n+i, j)); 79 } 80 } 81 } 82 83 // Check that the matrix is symmetric under time-reversal 84 bool is_t_symmetric(const double thresh = 1.0e-8) const { 85 const double err = check_t_symmetry(); 86 return err < thresh; 87 } 88 check_t_symmetry()89 double check_t_symmetry() const { 90 assert(mdim()%2 == 0 && ndim()%2 == 0); 91 const int m = mdim()/2; 92 const int n = ndim()/2; 93 94 std::shared_ptr<ZMatrix> u = get_submatrix(0, 0, n, m); 95 std::shared_ptr<ZMatrix> v = get_submatrix(n, 0, n, m); 96 97 u->ax_plus_y(-1.0, *get_submatrix(n, m, n, m)->get_conjg()); 98 v->ax_plus_y( 1.0, *get_submatrix(0, m, n, m)->get_conjg()); 99 100 const double out = u->rms() + v->rms(); 101 return out; 102 } 103 104 }; 105 106 } 107 108 #endif 109