Line data Source code
1 : /* +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ 2 : Copyright (c) 2014-2017 The plumed team 3 : (see the PEOPLE file at the root of the distribution for a list of names) 4 : 5 : See http://www.plumed.org for more information. 6 : 7 : This file is part of plumed, version 2. 8 : 9 : plumed is free software: you can redistribute it and/or modify 10 : it under the terms of the GNU Lesser General Public License as published by 11 : the Free Software Foundation, either version 3 of the License, or 12 : (at your option) any later version. 13 : 14 : plumed is distributed in the hope that it will be useful, 15 : but WITHOUT ANY WARRANTY; without even the implied warranty of 16 : MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 17 : GNU Lesser General Public License for more details. 18 : 19 : You should have received a copy of the GNU Lesser General Public License 20 : along with plumed. If not, see <http://www.gnu.org/licenses/>. 21 : +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ */ 22 : #include "MatrixOperationBase.h" 23 : #include "core/ActionRegister.h" 24 : 25 : //+PLUMEDOC ANALYSIS DIAGONALIZE 26 : /* 27 : Calculate the eigenvalues and eigenvectors of a square matrix 28 : 29 : \par Examples 30 : 31 : */ 32 : //+ENDPLUMEDOC 33 : 34 : namespace PLMD { 35 : namespace matrixtools { 36 : 37 : class DiagonalizeMatrix : public MatrixOperationBase { 38 : private: 39 : std::vector<unsigned> desired_vectors; 40 : Matrix<double> mymatrix; 41 : std::vector<double> eigvals; 42 : Matrix<double> eigvecs; 43 : public: 44 : static void registerKeywords( Keywords& keys ); 45 : /// Constructor 46 : explicit DiagonalizeMatrix(const ActionOptions&); 47 : /// This is required to set the number of derivatives for the eigenvalues 48 276 : unsigned getNumberOfDerivatives() override { return getPntrToArgument(0)->getNumberOfValues(); } 49 : /// 50 : void prepare() override ; 51 : /// 52 : void calculate() override ; 53 : /// 54 : double getForceOnMatrixElement( const unsigned& jrow, const unsigned& krow ) const override; 55 : }; 56 : 57 : PLUMED_REGISTER_ACTION(DiagonalizeMatrix,"DIAGONALIZE") 58 : 59 41 : void DiagonalizeMatrix::registerKeywords( Keywords& keys ) { 60 41 : MatrixOperationBase::registerKeywords( keys ); 61 82 : keys.add("compulsory","VECTORS","all","the eigenvalues and vectors that you would like to calculate. 1=largest, 2=second largest and so on"); 62 82 : keys.addOutputComponent("vals","default","the eigevalues of the input matrix"); 63 82 : keys.addOutputComponent("vecs","default","the eigenvectors of the input matrix"); 64 41 : } 65 : 66 22 : DiagonalizeMatrix::DiagonalizeMatrix(const ActionOptions& ao): 67 : Action(ao), 68 22 : MatrixOperationBase(ao) 69 : { 70 22 : if( getPntrToArgument(0)->getShape()[0]!=getPntrToArgument(0)->getShape()[1] ) error("input matrix should be square"); 71 : 72 44 : std::vector<std::string> eigv; parseVector("VECTORS",eigv); 73 22 : if( eigv.size()>1 ) { 74 10 : Tools::interpretRanges(eigv); desired_vectors.resize( eigv.size() ); 75 30 : for(unsigned i=0; i<eigv.size(); ++i) Tools::convert( eigv[i], desired_vectors[i] ); 76 : } else { 77 12 : if( eigv.size()==0 ) error("missing input to VECTORS keyword"); 78 : unsigned ivec; 79 12 : if( eigv[0]=="all" ) { 80 7 : desired_vectors.resize( getPntrToArgument(0)->getShape()[0] ); 81 21 : for(unsigned i=0; i<desired_vectors.size(); ++i) desired_vectors[i] = i + 1; 82 : } else { 83 5 : Tools::convert( eigv[0], ivec ); 84 5 : desired_vectors.resize(1); desired_vectors[0]=ivec; 85 : } 86 : } 87 : 88 22 : std::string num; std::vector<unsigned> eval_shape(0); 89 22 : std::vector<unsigned> evec_shape(1); evec_shape[0] = getPntrToArgument(0)->getShape()[0]; 90 61 : for(unsigned i=0; i<desired_vectors.size(); ++i) { 91 39 : Tools::convert( desired_vectors[i], num ); 92 78 : addComponent( "vals-" + num, eval_shape ); componentIsNotPeriodic( "vals-" + num ); 93 78 : addComponent( "vecs-" + num, evec_shape ); componentIsNotPeriodic( "vecs-" + num ); 94 : // Make sure eigenvalues are always stored 95 39 : getPntrToComponent( 2*i+1 )->buildDataStore(); 96 : } 97 : 98 22 : std::vector<unsigned> eigvecs_shape(2); eigvecs_shape[0]=eigvecs_shape[1]=getPntrToArgument(0)->getShape()[0]; 99 22 : mymatrix.resize( eigvecs_shape[0], eigvecs_shape[1] ); eigvals.resize( eigvecs_shape[0] ); eigvecs.resize( eigvecs_shape[0], eigvecs_shape[1] ); 100 22 : } 101 : 102 124 : void DiagonalizeMatrix::prepare() { 103 124 : std::vector<unsigned> shape(1); shape[0]=getPntrToArgument(0)->getShape()[0]; 104 303 : for(unsigned i=0; i<desired_vectors.size(); ++i) { 105 179 : if( getPntrToComponent( 2*i+1 )->getShape()[0]!=shape[0] ) getPntrToComponent( 2*i+1 )->setShape( shape ); 106 : } 107 : 108 124 : } 109 : 110 122 : void DiagonalizeMatrix::calculate() { 111 122 : if( getPntrToArgument(0)->getShape()[0]==0 ) return ; 112 : // Resize stuff that might need resizing 113 : unsigned nvals=getPntrToArgument(0)->getShape()[0]; 114 127 : if( eigvals.size()!=nvals ) { mymatrix.resize( nvals, nvals ); eigvals.resize( nvals ); eigvecs.resize( nvals, nvals ); } 115 : 116 : // Retrieve the matrix from input 117 122 : retrieveFullMatrix( mymatrix ); 118 : // Now diagonalize the matrix 119 122 : diagMat( mymatrix, eigvals, eigvecs ); 120 : // And set the eigenvalues and eigenvectors 121 297 : for(unsigned i=0; i<desired_vectors.size(); ++i) { 122 175 : getPntrToComponent(2*i)->set( eigvals[ mymatrix.ncols()-desired_vectors[i]] ); 123 175 : Value* evec_out = getPntrToComponent(2*i+1); unsigned vreq = mymatrix.ncols()-desired_vectors[i]; 124 4526 : for(unsigned j=0; j<mymatrix.ncols(); ++j) evec_out->set( j, eigvecs( vreq, j ) ); 125 : } 126 : } 127 : 128 12495 : double DiagonalizeMatrix::getForceOnMatrixElement( const unsigned& jrow, const unsigned& kcol ) const { 129 : double ff = 0; 130 31654 : for(unsigned i=0; i<desired_vectors.size(); ++i) { 131 : // Deal with forces on eigenvalues 132 19159 : if( getConstPntrToComponent(2*i)->forcesWereAdded() ) { 133 8330 : unsigned ncol = mymatrix.ncols()-desired_vectors[i]; 134 8330 : ff += getConstPntrToComponent(2*i)->getForce(0)*eigvecs(ncol,jrow)*eigvecs(ncol,kcol); 135 : } 136 : // And forces on eigenvectors 137 19159 : if( !getConstPntrToComponent(2*i+1)->forcesWereAdded() ) continue; 138 : 139 7497 : unsigned ncol = mymatrix.ncols()-desired_vectors[i]; 140 106624 : for(unsigned n=0; n<mymatrix.nrows(); ++n) { 141 : double tmp2 = 0; 142 1446088 : for(unsigned m=0; m<mymatrix.nrows(); ++m) { 143 1346961 : if( m==ncol ) continue; 144 1247834 : tmp2 += eigvecs(m,n)*eigvecs(m,jrow)*eigvecs(ncol,kcol) / (eigvals[ncol]-eigvals[m]); 145 : } 146 99127 : ff += getConstPntrToComponent(2*i+1)->getForce(n) * tmp2; 147 : } 148 : } 149 12495 : return ff; 150 : } 151 : 152 : } 153 : }