Line data Source code
1 : /* +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
2 : Copyright (c) 2013-2019 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 6459 : 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 21 : MultiDomainRMSD::~MultiDomainRMSD() {
39 56 : for(unsigned i=0; i<domains.size(); ++i) delete domains[i];
40 14 : }
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 : std::vector<Vector> positions; std::vector<double> align, displace;
47 14 : std::string num; blocks.resize( nblocks+1 ); blocks[0]=0;
48 49 : 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 14 : parse("LOWER_CUTOFF",lower,true);
52 14 : parse("UPPER_CUTOFF",upper,true);
53 14 : bool nopbc=false; parseFlag("NOPBC",nopbc);
54 :
55 35 : for(unsigned i=1; i<=nblocks; ++i) {
56 14 : Tools::convert(i,num);
57 28 : 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 28 : domains.push_back( metricRegister().create<SingleDomainRMSD>( ftype ) );
64 42 : positions.resize( blocks[i] - blocks[i-1] );
65 28 : align.resize( blocks[i] - blocks[i-1] );
66 28 : displace.resize( blocks[i] - blocks[i-1] );
67 : unsigned n=0;
68 174 : for(unsigned j=blocks[i-1]; j<blocks[i]; ++j) {
69 219 : positions[n]=pdb.getPositions()[j];
70 146 : align[n]=pdb.getOccupancy()[j];
71 146 : 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 28 : double ww=0; parse("WEIGHT"+num, ww, true );
79 28 : 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 7 : 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 111 : 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 222 : for(unsigned i=0; i<domains.size(); ++i) {
95 : // Must extract appropriate positions here
96 222 : mypos.resize( blocks[i+1] - blocks[i] );
97 118 : if( myder.calcUsingPCAOption() ) domains[i]->setupPCAStorage( tder );
98 1285 : unsigned n=0; for(unsigned j=blocks[i]; j<blocks[i+1]; ++j) { tder.setAtomIndex(n,j); mypos[n]=pos[j]; n++; }
99 1211 : for(unsigned k=n; k<getNumberOfAtoms(); ++k) tder.setAtomIndex(k,3*pos.size()+10);
100 : // This actually does the calculation
101 148 : 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 : unsigned n=0;
107 66 : if( tder.centeredpos.size()>0 ) myder.rot[i]=tder.rot[0];
108 396 : for(unsigned j=blocks[i]; j<blocks[i+1]; ++j) {
109 462 : 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 1001 : 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 37 : 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 : bool enabled=true;
138 16 : 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 : mypack.switchOnPCAOption();
147 4 : mypack.displacement.resize( getNumberOfAtoms() );
148 4 : mypack.centeredpos.resize( getNumberOfAtoms() );
149 4 : mypack.DRotDPos.resize(3,3); mypack.rot.resize( domains.size() );
150 26 : 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 0 : void MultiDomainRMSD::extractAtomicDisplacement( const std::vector<Vector>& pos, std::vector<Vector>& direction ) const {
164 : std::vector<Vector> mypos, mydir;
165 0 : for(unsigned i=0; i<domains.size(); ++i) {
166 : // Must extract appropriate positions here
167 0 : mypos.resize( blocks[i+1] - blocks[i] ); mydir.resize( blocks[i+1] - blocks[i] );
168 0 : unsigned n=0; for(unsigned j=blocks[i]; j<blocks[i+1]; ++j) { mypos[n]=pos[j]; n++; }
169 : // Do the calculation
170 0 : domains[i]->extractAtomicDisplacement( mypos, mydir );
171 : // Extract the direction
172 0 : n=0; for(unsigned j=blocks[i]; j<blocks[i+1]; ++j) { direction[j]=weights[i]*mydir[n]; n++; }
173 : }
174 0 : }
175 :
176 44 : double MultiDomainRMSD::projectAtomicDisplacementOnVector( const bool& normalized, const std::vector<Vector>& vecs, ReferenceValuePack& mypack ) const {
177 44 : double totd=0.; std::vector<Vector> tvecs; mypack.clear();
178 132 : MultiValue tvals( 1, mypack.getNumberOfDerivatives() ); ReferenceValuePack tder( 0, getNumberOfAtoms(), tvals );
179 352 : for(unsigned i=0; i<domains.size(); ++i) {
180 : // Must extract appropriate positions here
181 352 : tvecs.resize( blocks[i+1] - blocks[i] ); domains[i]->setupPCAStorage( tder );
182 88 : if( tder.centeredpos.size()>0 ) {
183 572 : for(unsigned p=0; p<3; ++p) for(unsigned q=0; q<3; ++q) tder.DRotDPos(p,q).resize( tvecs.size() );
184 : }
185 : // Extract information from storage pack and put in local pack
186 132 : if( tder.centeredpos.size()>0 ) tder.rot[0]=mypack.rot[i];
187 : unsigned n=0;
188 792 : for(unsigned j=blocks[i]; j<blocks[i+1]; ++j) {
189 1232 : tder.setAtomIndex(n,j); tvecs[n] = vecs[j]; tder.displacement[n]=mypack.displacement[j] / weights[i];
190 308 : if( tder.centeredpos.size()>0 ) {
191 154 : tder.centeredpos[n]=mypack.centeredpos[j];
192 2002 : for(unsigned p=0; p<3; ++p) for(unsigned q=0; q<3; ++q) tder.DRotDPos(p,q)[n]=mypack.DRotDPos(p,q)[j];
193 : }
194 308 : n++;
195 : }
196 1012 : for(unsigned k=n; k<getNumberOfAtoms(); ++k) tder.setAtomIndex(k,3*vecs.size()+10);
197 :
198 : // Do the calculations
199 176 : totd += weights[i]*domains[i]->projectAtomicDisplacementOnVector( normalized, tvecs, tder );
200 :
201 : // And derivatives
202 88 : mypack.copyScaledDerivatives( 0, weights[i], tvals );
203 : }
204 44 : if( !mypack.updateComplete() ) mypack.updateDynamicLists();
205 :
206 44 : return totd;
207 : }
208 :
209 4839 : }
|
