Line data Source code
1 : /* +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ 2 : Copyright (c) 2013-2023 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 "SingleDomainRMSD.h" 23 : #include "tools/PDB.h" 24 : #include "DRMSD.h" 25 : 26 : namespace PLMD { 27 : 28 494 : SingleDomainRMSD::SingleDomainRMSD( const ReferenceConfigurationOptions& ro ): 29 : ReferenceConfiguration(ro), 30 494 : ReferenceAtoms(ro) 31 : { 32 494 : } 33 : 34 452 : void SingleDomainRMSD::readReference( const PDB& pdb ) { 35 452 : readAtomsFromPDB( pdb ); 36 : double wa=0, wd=0; 37 8327 : for(unsigned i=0; i<pdb.size(); ++i) { wa+=align[i]; wd+=displace[i]; } 38 : 39 452 : if(wa>epsilon) { 40 450 : double w=1.0/wa; 41 8295 : for(unsigned i=0; i<pdb.size(); ++i) align[i] *= w; 42 : } else { 43 2 : double w=1.0/pdb.size(); 44 32 : for(unsigned i=0; i<pdb.size(); ++i) align[i] = w; 45 : } 46 : 47 452 : if(wd>epsilon) { 48 450 : double w=1.0/wd; 49 8295 : for(unsigned i=0; i<pdb.size(); ++i) displace[i] *= w; 50 : } else { 51 2 : double w=1.0/pdb.size(); 52 32 : for(unsigned i=0; i<pdb.size(); ++i) displace[i] = w; 53 : } 54 : 55 452 : Vector center; 56 8327 : for(unsigned i=0; i<pdb.size(); ++i) { 57 7875 : center+=reference_atoms[i]*align[i]; 58 : } 59 8327 : for(unsigned i=0; i<pdb.size(); ++i) reference_atoms[i]-=center; 60 452 : } 61 : 62 45 : void SingleDomainRMSD::setReferenceAtoms( const std::vector<Vector>& conf, const std::vector<double>& align_in, const std::vector<double>& displace_in ) { 63 45 : reference_atoms.resize( conf.size() ); align.resize( conf.size() ); 64 45 : displace.resize( conf.size() ); atom_der_index.resize( conf.size() ); 65 : double wa=0, wd=0; 66 1154 : for(unsigned i=0; i<conf.size(); ++i) { wa+=align_in[i]; wd+=displace_in[i]; } 67 : 68 45 : if(wa>epsilon) { 69 44 : double w=1.0/wa; 70 1145 : for(unsigned i=0; i<conf.size(); ++i) align[i] = align_in[i] * w; 71 : } else { 72 1 : double w=1.0/conf.size(); 73 9 : for(unsigned i=0; i<conf.size(); ++i) align[i] = w; 74 : } 75 : 76 45 : if(wd>epsilon) { 77 44 : double w=1.0/wd; 78 1146 : for(unsigned i=0; i<conf.size(); ++i) displace[i] = displace_in[i] * w; 79 : } else { 80 1 : double w=1.0/conf.size(); 81 8 : for(unsigned i=0; i<conf.size(); ++i) displace[i] = w; 82 : } 83 : 84 45 : Vector center; 85 1154 : for(unsigned i=0; i<conf.size(); ++i) { 86 1109 : center+=conf[i]*align[i]; atom_der_index[i]=i; 87 : } 88 1154 : for(unsigned i=0; i<conf.size(); ++i) reference_atoms[i]=conf[i]-center; 89 45 : setupRMSDObject(); 90 45 : } 91 : 92 58046 : double SingleDomainRMSD::calculate( const std::vector<Vector>& pos, const Pbc& pbc, ReferenceValuePack& myder, const bool& squared ) const { 93 58046 : return calc( pos, pbc, myder, squared ); 94 : } 95 : 96 140941 : double SingleDomainRMSD::calc( const std::vector<Vector>& pos, const Pbc& pbc, const std::vector<Value*>& vals, const std::vector<double>& arg, 97 : ReferenceValuePack& myder, const bool& squared ) const { 98 : plumed_dbg_assert( vals.size()==0 && pos.size()==getNumberOfAtoms() && arg.size()==0 ); 99 140941 : return calc( pos, pbc, myder, squared ); 100 : } 101 : 102 : }