Line data Source code
1 : /* +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
2 : Copyright (c) 2013-2018 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 "MultiDomainRMSD.h"
23 : #include "SingleDomainRMSD.h"
24 : #include "MetricRegister.h"
25 : #include "tools/PDB.h"
26 :
27 : namespace PLMD {
28 :
29 2530 : PLUMED_REGISTER_METRIC(MultiDomainRMSD,"MULTI")
30 :
31 7 : MultiDomainRMSD::MultiDomainRMSD( const ReferenceConfigurationOptions& ro ):
32 : ReferenceConfiguration(ro),
33 : ReferenceAtoms(ro),
34 7 : ftype(ro.getMultiRMSDType())
35 : {
36 7 : }
37 :
38 28 : MultiDomainRMSD::~MultiDomainRMSD() {
39 7 : for(unsigned i=0; i<domains.size(); ++i) delete domains[i];
40 21 : }
41 :
42 7 : void MultiDomainRMSD::read( const PDB& pdb ) {
43 7 : unsigned nblocks = pdb.getNumberOfAtomBlocks();
44 7 : if( nblocks<2 ) error("multidomain RMSD only has one block of atoms");
45 :
46 14 : std::vector<Vector> positions; std::vector<double> align, displace;
47 14 : std::string num; blocks.resize( nblocks+1 ); blocks[0]=0;
48 7 : for(unsigned i=0; i<nblocks; ++i) blocks[i+1]=pdb.getAtomBlockEnds()[i];
49 :
50 7 : double lower=0.0, upper=std::numeric_limits<double>::max( );
51 7 : parse("LOWER_CUTOFF",lower,true);
52 7 : parse("UPPER_CUTOFF",upper,true);
53 7 : bool nopbc=false; parseFlag("NOPBC",nopbc);
54 :
55 21 : for(unsigned i=1; i<=nblocks; ++i) {
56 14 : Tools::convert(i,num);
57 14 : if( ftype=="RMSD" ) {
58 0 : parse("TYPE"+num, ftype );
59 0 : parse("LOWER_CUTOFF"+num,lower,true);
60 0 : parse("UPPER_CUTOFF"+num,upper,true);
61 0 : nopbc=false; parseFlag("NOPBC"+num,nopbc);
62 : }
63 14 : domains.push_back( metricRegister().create<SingleDomainRMSD>( ftype ) );
64 14 : positions.resize( blocks[i] - blocks[i-1] );
65 14 : align.resize( blocks[i] - blocks[i-1] );
66 14 : displace.resize( blocks[i] - blocks[i-1] );
67 14 : unsigned n=0;
68 87 : for(unsigned j=blocks[i-1]; j<blocks[i]; ++j) {
69 73 : positions[n]=pdb.getPositions()[j];
70 73 : align[n]=pdb.getOccupancy()[j];
71 73 : displace[n]=pdb.getBeta()[j];
72 73 : n++;
73 : }
74 14 : domains[i-1]->setBoundsOnDistances( !nopbc, lower, upper );
75 14 : domains[i-1]->setReferenceAtoms( positions, align, displace );
76 14 : domains[i-1]->setupRMSDObject();
77 :
78 14 : double ww=0; parse("WEIGHT"+num, ww, true );
79 14 : if( ww==0 ) weights.push_back( 1.0 );
80 0 : else weights.push_back( ww );
81 : }
82 : // And set the atom numbers for this object
83 14 : setAtomNumbers( pdb.getAtomNumbers() );
84 7 : }
85 :
86 0 : void MultiDomainRMSD::setReferenceAtoms( const std::vector<Vector>& conf, const std::vector<double>& align_in, const std::vector<double>& displace_in ) {
87 0 : plumed_error();
88 : }
89 :
90 37 : double MultiDomainRMSD::calculate( const std::vector<Vector>& pos, const Pbc& pbc, ReferenceValuePack& myder, const bool& squared ) const {
91 74 : double totd=0.; Tensor tvirial; std::vector<Vector> mypos; MultiValue tvals( 1, 3*pos.size()+9 );
92 74 : ReferenceValuePack tder( 0, getNumberOfAtoms(), tvals ); myder.clear();
93 :
94 111 : for(unsigned i=0; i<domains.size(); ++i) {
95 : // Must extract appropriate positions here
96 74 : mypos.resize( blocks[i+1] - blocks[i] );
97 74 : if( myder.calcUsingPCAOption() ) domains[i]->setupPCAStorage( tder );
98 74 : unsigned n=0; for(unsigned j=blocks[i]; j<blocks[i+1]; ++j) { tder.setAtomIndex(n,j); mypos[n]=pos[j]; n++; }
99 74 : for(unsigned k=n; k<getNumberOfAtoms(); ++k) tder.setAtomIndex(k,3*pos.size()+10);
100 : // This actually does the calculation
101 74 : totd += weights[i]*domains[i]->calculate( mypos, pbc, tder, true );
102 : // Now merge the derivative
103 74 : myder.copyScaledDerivatives( 0, weights[i], tvals );
104 : // If PCA copy PCA stuff
105 74 : if( myder.calcUsingPCAOption() ) {
106 44 : unsigned n=0;
107 44 : if( tder.centeredpos.size()>0 ) myder.rot[i]=tder.rot[0];
108 198 : for(unsigned j=blocks[i]; j<blocks[i+1]; ++j) {
109 154 : myder.displacement[j]=weights[i]*tder.displacement[n]; // Multiplication by weights here ensures that normalisation is done correctly
110 154 : if( tder.centeredpos.size()>0 ) {
111 77 : myder.centeredpos[j]=tder.centeredpos[n];
112 77 : for(unsigned p=0; p<3; ++p) for(unsigned q=0; q<3; ++q) myder.DRotDPos(p,q)[j]=tder.DRotDPos(p,q)[n];
113 : }
114 154 : n++;
115 : }
116 : }
117 : // Make sure virial status is set correctly in output derivative pack
118 : // This is only done here so I do this by using class friendship
119 74 : if( tder.virialWasSet() ) myder.boxWasSet=true;
120 : }
121 37 : if( !myder.updateComplete() ) myder.updateDynamicLists();
122 :
123 37 : if( !squared ) {
124 15 : totd=sqrt(totd); double xx=0.5/totd;
125 15 : myder.scaleAllDerivatives( xx );
126 : }
127 74 : return totd;
128 : }
129 :
130 22 : double MultiDomainRMSD::calc( const std::vector<Vector>& pos, const Pbc& pbc, const std::vector<Value*>& vals, const std::vector<double>& arg,
131 : ReferenceValuePack& myder, const bool& squared ) const {
132 : plumed_dbg_assert( vals.size()==0 && pos.size()==getNumberOfAtoms() && arg.size()==0 );
133 22 : return calculate( pos, pbc, myder, squared );
134 : }
135 :
136 2 : bool MultiDomainRMSD::pcaIsEnabledForThisReference() {
137 2 : bool enabled=true;
138 6 : for(unsigned i=0; i<domains.size(); ++i) {
139 4 : if( !domains[i]->pcaIsEnabledForThisReference() ) enabled=false;
140 : }
141 2 : return enabled;
142 : }
143 :
144 2 : void MultiDomainRMSD::setupPCAStorage( ReferenceValuePack& mypack ) {
145 : plumed_dbg_assert( pcaIsEnabledForThisReference() );
146 2 : mypack.switchOnPCAOption();
147 2 : mypack.displacement.resize( getNumberOfAtoms() );
148 2 : mypack.centeredpos.resize( getNumberOfAtoms() );
149 2 : mypack.DRotDPos.resize(3,3); mypack.rot.resize( domains.size() );
150 2 : for(unsigned i=0; i<3; ++i) for(unsigned j=0; j<3; ++j) mypack.DRotDPos(i,j).resize( getNumberOfAtoms() );
151 2 : }
152 :
153 : // Vector MultiDomainRMSD::getAtomicDisplacement( const unsigned& iatom ){
154 : // for(unsigned i=0;i<domains.size();++i){
155 : // unsigned n=0;
156 : // for(unsigned j=blocks[i];j<blocks[i+1];++j){
157 : // if( j==iatom ) return weights[i]*domains[i]->getAtomicDisplacement(n);
158 : // n++;
159 : // }
160 : // }
161 : // }
162 :
163 44 : double MultiDomainRMSD::projectAtomicDisplacementOnVector( const unsigned& iv, const Matrix<Vector>& vecs, const std::vector<Vector>& pos, ReferenceValuePack& mypack ) const {
164 88 : double totd=0.; Matrix<Vector> tvecs; std::vector<Vector> mypos; mypack.clear();
165 88 : MultiValue tvals( 1, mypack.getNumberOfDerivatives() ); ReferenceValuePack tder( 0, getNumberOfAtoms(), tvals );
166 132 : for(unsigned i=0; i<domains.size(); ++i) {
167 : // Must extract appropriate positions here
168 88 : mypos.resize( blocks[i+1] - blocks[i] ); tvecs.resize( 1, blocks[i+1] - blocks[i] );
169 88 : domains[i]->setupPCAStorage( tder );
170 88 : if( tder.centeredpos.size()>0 ) {
171 44 : for(unsigned p=0; p<3; ++p) for(unsigned q=0; q<3; ++q) tder.DRotDPos(p,q).resize( mypos.size() );
172 : }
173 : // Extract information from storage pack and put in local pack
174 88 : if( tder.centeredpos.size()>0 ) tder.rot[0]=mypack.rot[i];
175 88 : unsigned n=0;
176 396 : for(unsigned j=blocks[i]; j<blocks[i+1]; ++j) {
177 308 : mypos[n]=pos[j]; tder.setAtomIndex(n,j); tvecs( 0, n ) = vecs( iv, j );
178 308 : tder.displacement[n]=mypack.displacement[j] / weights[i];
179 308 : if( tder.centeredpos.size()>0 ) {
180 154 : tder.centeredpos[n]=mypack.centeredpos[j];
181 154 : for(unsigned p=0; p<3; ++p) for(unsigned q=0; q<3; ++q) tder.DRotDPos(p,q)[n]=mypack.DRotDPos(p,q)[j];
182 : }
183 308 : n++;
184 : }
185 88 : for(unsigned k=n; k<getNumberOfAtoms(); ++k) tder.setAtomIndex(k,3*pos.size()+10);
186 :
187 : // Do the calculations
188 88 : totd += weights[i]*domains[i]->projectAtomicDisplacementOnVector( 0, tvecs, mypos, tder );
189 :
190 : // And derivatives
191 88 : mypack.copyScaledDerivatives( 0, weights[i], tvals );
192 : }
193 44 : if( !mypack.updateComplete() ) mypack.updateDynamicLists();
194 :
195 88 : return totd;
196 : }
197 :
198 2523 : }
|
