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 "MultiDomainRMSD.h"
23 : #include "SingleDomainRMSD.h"
24 : #include "MetricRegister.h"
25 : #include "tools/PDB.h"
26 :
27 : namespace PLMD {
28 :
29 13795 : PLUMED_REGISTER_METRIC(MultiDomainRMSD,"MULTI")
30 :
31 5 : MultiDomainRMSD::MultiDomainRMSD( const ReferenceConfigurationOptions& ro ):
32 : ReferenceConfiguration(ro),
33 : ReferenceAtoms(ro),
34 5 : ftype(ro.getMultiRMSDType()) {
35 5 : }
36 :
37 5 : void MultiDomainRMSD::read( const PDB& pdb ) {
38 5 : unsigned nblocks = pdb.getNumberOfAtomBlocks();
39 5 : if( nblocks<2 ) {
40 0 : error("multidomain RMSD only has one block of atoms");
41 : }
42 :
43 : std::vector<Vector> positions;
44 : std::vector<double> align, displace;
45 : std::string num;
46 5 : blocks.resize( nblocks+1 );
47 5 : blocks[0]=0;
48 15 : for(unsigned i=0; i<nblocks; ++i) {
49 10 : blocks[i+1]=pdb.getAtomBlockEnds()[i];
50 : }
51 :
52 : double tmp, lower=0.0, upper=std::numeric_limits<double>::max( );
53 10 : if( pdb.getArgumentValue("LOWER_CUTOFF",tmp) ) {
54 0 : lower=tmp;
55 : }
56 10 : if( pdb.getArgumentValue("UPPER_CUTOFF",tmp) ) {
57 0 : upper=tmp;
58 : }
59 5 : bool nopbc=pdb.hasFlag("NOPBC");
60 :
61 5 : domains.resize(0);
62 5 : weights.resize(0);
63 15 : for(unsigned i=1; i<=nblocks; ++i) {
64 10 : Tools::convert(i,num);
65 10 : if( ftype=="RMSD" ) {
66 : // parse("TYPE"+num, ftype );
67 : lower=0.0;
68 : upper=std::numeric_limits<double>::max( );
69 0 : if( pdb.getArgumentValue("LOWER_CUTOFF"+num,tmp) ) {
70 0 : lower=tmp;
71 : }
72 0 : if( pdb.getArgumentValue("UPPER_CUTOFF"+num,tmp) ) {
73 0 : upper=tmp;
74 : }
75 0 : nopbc=pdb.hasFlag("NOPBC");
76 : }
77 10 : domains.emplace_back( metricRegister().create<SingleDomainRMSD>( ftype ) );
78 10 : positions.resize( blocks[i] - blocks[i-1] );
79 10 : align.resize( blocks[i] - blocks[i-1] );
80 10 : displace.resize( blocks[i] - blocks[i-1] );
81 : unsigned n=0;
82 69 : for(unsigned j=blocks[i-1]; j<blocks[i]; ++j) {
83 59 : positions[n]=pdb.getPositions()[j];
84 59 : align[n]=pdb.getOccupancy()[j];
85 59 : displace[n]=pdb.getBeta()[j];
86 59 : n++;
87 : }
88 10 : domains[i-1]->setBoundsOnDistances( !nopbc, lower, upper );
89 10 : domains[i-1]->setReferenceAtoms( positions, align, displace );
90 10 : domains[i-1]->setupRMSDObject();
91 :
92 10 : double ww=0;
93 20 : if( !pdb.getArgumentValue("WEIGHT"+num,ww) ) {
94 10 : weights.push_back( 1.0 );
95 : } else {
96 0 : weights.push_back( ww );
97 : }
98 : }
99 : // And set the atom numbers for this object
100 5 : indices.resize(0);
101 5 : atom_der_index.resize(0);
102 64 : for(unsigned i=0; i<pdb.size(); ++i) {
103 59 : indices.push_back( pdb.getAtomNumbers()[i] );
104 59 : atom_der_index.push_back(i);
105 : }
106 : // setAtomNumbers( pdb.getAtomNumbers() );
107 5 : }
108 :
109 0 : void MultiDomainRMSD::setReferenceAtoms( const std::vector<Vector>& conf, const std::vector<double>& align_in, const std::vector<double>& displace_in ) {
110 0 : plumed_error();
111 : }
112 :
113 37 : double MultiDomainRMSD::calculate( const std::vector<Vector>& pos, const Pbc& pbc, ReferenceValuePack& myder, const bool& squared ) const {
114 : double totd=0.;
115 37 : Tensor tvirial;
116 : std::vector<Vector> mypos;
117 37 : MultiValue tvals( 1, 3*pos.size()+9 );
118 37 : ReferenceValuePack tder( 0, getNumberOfAtoms(), tvals );
119 37 : myder.clear();
120 :
121 111 : for(unsigned i=0; i<domains.size(); ++i) {
122 : // Must extract appropriate positions here
123 74 : mypos.resize( blocks[i+1] - blocks[i] );
124 74 : if( myder.calcUsingPCAOption() ) {
125 44 : domains[i]->setupPCAStorage( tder );
126 : }
127 : unsigned n=0;
128 453 : for(unsigned j=blocks[i]; j<blocks[i+1]; ++j) {
129 : tder.setAtomIndex(n,j);
130 379 : mypos[n]=pos[j];
131 379 : n++;
132 : }
133 379 : for(unsigned k=n; k<getNumberOfAtoms(); ++k) {
134 379 : tder.setAtomIndex(k,3*pos.size()+10);
135 : }
136 : // This actually does the calculation
137 74 : totd += weights[i]*domains[i]->calculate( mypos, pbc, tder, true );
138 : // Now merge the derivative
139 74 : myder.copyScaledDerivatives( 0, weights[i], tvals );
140 : // If PCA copy PCA stuff
141 74 : if( myder.calcUsingPCAOption() ) {
142 : unsigned n=0;
143 44 : if( tder.centeredpos.size()>0 ) {
144 22 : myder.rot[i]=tder.rot[0];
145 : }
146 198 : for(unsigned j=blocks[i]; j<blocks[i+1]; ++j) {
147 154 : myder.displacement[j]=weights[i]*tder.displacement[n]; // Multiplication by weights here ensures that normalisation is done correctly
148 154 : if( tder.centeredpos.size()>0 ) {
149 77 : myder.centeredpos[j]=tder.centeredpos[n];
150 308 : for(unsigned p=0; p<3; ++p)
151 924 : for(unsigned q=0; q<3; ++q) {
152 693 : myder.DRotDPos(p,q)[j]=tder.DRotDPos(p,q)[n];
153 : }
154 : }
155 154 : n++;
156 : }
157 : }
158 : // Make sure virial status is set correctly in output derivative pack
159 : // This is only done here so I do this by using class friendship
160 74 : if( tder.virialWasSet() ) {
161 10 : myder.boxWasSet=true;
162 : }
163 : }
164 37 : if( !myder.updateComplete() ) {
165 37 : myder.updateDynamicLists();
166 : }
167 :
168 37 : if( !squared ) {
169 15 : totd=std::sqrt(totd);
170 15 : double xx=0.5/totd;
171 15 : myder.scaleAllDerivatives( xx );
172 : }
173 37 : return totd;
174 37 : }
175 :
176 22 : double MultiDomainRMSD::calc( const std::vector<Vector>& pos, const Pbc& pbc, const std::vector<Value*>& vals, const std::vector<double>& arg,
177 : ReferenceValuePack& myder, const bool& squared ) const {
178 : plumed_dbg_assert( vals.size()==0 && pos.size()==getNumberOfAtoms() && arg.size()==0 );
179 22 : return calculate( pos, pbc, myder, squared );
180 : }
181 :
182 2 : bool MultiDomainRMSD::pcaIsEnabledForThisReference() {
183 : bool enabled=true;
184 6 : for(unsigned i=0; i<domains.size(); ++i) {
185 4 : if( !domains[i]->pcaIsEnabledForThisReference() ) {
186 : enabled=false;
187 : }
188 : }
189 2 : return enabled;
190 : }
191 :
192 2 : void MultiDomainRMSD::setupPCAStorage( ReferenceValuePack& mypack ) {
193 : plumed_dbg_assert( pcaIsEnabledForThisReference() );
194 : mypack.switchOnPCAOption();
195 2 : mypack.displacement.resize( getNumberOfAtoms() );
196 2 : mypack.centeredpos.resize( getNumberOfAtoms() );
197 : mypack.DRotDPos.resize(3,3);
198 2 : mypack.rot.resize( domains.size() );
199 8 : for(unsigned i=0; i<3; ++i)
200 24 : for(unsigned j=0; j<3; ++j) {
201 18 : mypack.DRotDPos(i,j).resize( getNumberOfAtoms() );
202 : }
203 2 : }
204 :
205 : // Vector MultiDomainRMSD::getAtomicDisplacement( const unsigned& iatom ){
206 : // for(unsigned i=0;i<domains.size();++i){
207 : // unsigned n=0;
208 : // for(unsigned j=blocks[i];j<blocks[i+1];++j){
209 : // if( j==iatom ) return weights[i]*domains[i]->getAtomicDisplacement(n);
210 : // n++;
211 : // }
212 : // }
213 : // }
214 :
215 0 : void MultiDomainRMSD::extractAtomicDisplacement( const std::vector<Vector>& pos, std::vector<Vector>& direction ) const {
216 : std::vector<Vector> mypos, mydir;
217 0 : for(unsigned i=0; i<domains.size(); ++i) {
218 : // Must extract appropriate positions here
219 0 : mypos.resize( blocks[i+1] - blocks[i] );
220 0 : mydir.resize( blocks[i+1] - blocks[i] );
221 : unsigned n=0;
222 0 : for(unsigned j=blocks[i]; j<blocks[i+1]; ++j) {
223 0 : mypos[n]=pos[j];
224 0 : n++;
225 : }
226 : // Do the calculation
227 0 : domains[i]->extractAtomicDisplacement( mypos, mydir );
228 : // Extract the direction
229 : n=0;
230 0 : for(unsigned j=blocks[i]; j<blocks[i+1]; ++j) {
231 0 : direction[j]=weights[i]*mydir[n];
232 0 : n++;
233 : }
234 : }
235 0 : }
236 :
237 44 : double MultiDomainRMSD::projectAtomicDisplacementOnVector( const bool& normalized, const std::vector<Vector>& vecs, ReferenceValuePack& mypack ) const {
238 : double totd=0.;
239 : std::vector<Vector> tvecs;
240 44 : mypack.clear();
241 44 : MultiValue tvals( 1, mypack.getNumberOfDerivatives() );
242 44 : ReferenceValuePack tder( 0, getNumberOfAtoms(), tvals );
243 132 : for(unsigned i=0; i<domains.size(); ++i) {
244 : // Must extract appropriate positions here
245 88 : tvecs.resize( blocks[i+1] - blocks[i] );
246 88 : domains[i]->setupPCAStorage( tder );
247 88 : if( tder.centeredpos.size()>0 ) {
248 176 : for(unsigned p=0; p<3; ++p)
249 528 : for(unsigned q=0; q<3; ++q) {
250 396 : tder.DRotDPos(p,q).resize( tvecs.size() );
251 : }
252 : }
253 : // Extract information from storage pack and put in local pack
254 88 : if( tder.centeredpos.size()>0 ) {
255 44 : tder.rot[0]=mypack.rot[i];
256 : }
257 : unsigned n=0;
258 396 : for(unsigned j=blocks[i]; j<blocks[i+1]; ++j) {
259 : tder.setAtomIndex(n,j);
260 308 : tvecs[n] = vecs[j];
261 308 : tder.displacement[n]=mypack.displacement[j] / weights[i];
262 308 : if( tder.centeredpos.size()>0 ) {
263 154 : tder.centeredpos[n]=mypack.centeredpos[j];
264 616 : for(unsigned p=0; p<3; ++p)
265 1848 : for(unsigned q=0; q<3; ++q) {
266 1386 : tder.DRotDPos(p,q)[n]=mypack.DRotDPos(p,q)[j];
267 : }
268 : }
269 308 : n++;
270 : }
271 396 : for(unsigned k=n; k<getNumberOfAtoms(); ++k) {
272 308 : tder.setAtomIndex(k,3*vecs.size()+10);
273 : }
274 :
275 : // Do the calculations
276 88 : totd += weights[i]*domains[i]->projectAtomicDisplacementOnVector( normalized, tvecs, tder );
277 :
278 : // And derivatives
279 88 : mypack.copyScaledDerivatives( 0, weights[i], tvals );
280 : }
281 44 : if( !mypack.updateComplete() ) {
282 44 : mypack.updateDynamicLists();
283 : }
284 :
285 44 : return totd;
286 44 : }
287 :
288 : }
|
