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 10429 : 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 : {
36 5 : }
37 :
38 5 : void MultiDomainRMSD::read( const PDB& pdb ) {
39 5 : unsigned nblocks = pdb.getNumberOfAtomBlocks();
40 5 : if( nblocks<2 ) error("multidomain RMSD only has one block of atoms");
41 :
42 : std::vector<Vector> positions; std::vector<double> align, displace;
43 5 : std::string num; blocks.resize( nblocks+1 ); blocks[0]=0;
44 15 : for(unsigned i=0; i<nblocks; ++i) blocks[i+1]=pdb.getAtomBlockEnds()[i];
45 :
46 : double tmp, lower=0.0, upper=std::numeric_limits<double>::max( );
47 10 : if( pdb.getArgumentValue("LOWER_CUTOFF",tmp) ) lower=tmp;
48 10 : if( pdb.getArgumentValue("UPPER_CUTOFF",tmp) ) upper=tmp;
49 5 : bool nopbc=pdb.hasFlag("NOPBC");
50 :
51 5 : domains.resize(0); weights.resize(0);
52 15 : for(unsigned i=1; i<=nblocks; ++i) {
53 10 : Tools::convert(i,num);
54 10 : if( ftype=="RMSD" ) {
55 : // parse("TYPE"+num, ftype );
56 : lower=0.0; upper=std::numeric_limits<double>::max( );
57 0 : if( pdb.getArgumentValue("LOWER_CUTOFF"+num,tmp) ) lower=tmp;
58 0 : if( pdb.getArgumentValue("UPPER_CUTOFF"+num,tmp) ) upper=tmp;
59 0 : nopbc=pdb.hasFlag("NOPBC");
60 : }
61 10 : domains.emplace_back( metricRegister().create<SingleDomainRMSD>( ftype ) );
62 10 : positions.resize( blocks[i] - blocks[i-1] );
63 10 : align.resize( blocks[i] - blocks[i-1] );
64 10 : displace.resize( blocks[i] - blocks[i-1] );
65 : unsigned n=0;
66 69 : for(unsigned j=blocks[i-1]; j<blocks[i]; ++j) {
67 59 : positions[n]=pdb.getPositions()[j];
68 59 : align[n]=pdb.getOccupancy()[j];
69 59 : displace[n]=pdb.getBeta()[j];
70 59 : n++;
71 : }
72 10 : domains[i-1]->setBoundsOnDistances( !nopbc, lower, upper );
73 10 : domains[i-1]->setReferenceAtoms( positions, align, displace );
74 10 : domains[i-1]->setupRMSDObject();
75 :
76 10 : double ww=0;
77 20 : if( !pdb.getArgumentValue("WEIGHT"+num,ww) ) weights.push_back( 1.0 );
78 0 : else weights.push_back( ww );
79 : }
80 : // And set the atom numbers for this object
81 5 : indices.resize(0); atom_der_index.resize(0);
82 64 : for(unsigned i=0; i<pdb.size(); ++i) { indices.push_back( pdb.getAtomNumbers()[i] ); atom_der_index.push_back(i); }
83 : // setAtomNumbers( pdb.getAtomNumbers() );
84 5 : }
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 37 : double totd=0.; Tensor tvirial; std::vector<Vector> mypos; MultiValue tvals( 1, 3*pos.size()+9 );
92 37 : 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 453 : unsigned n=0; for(unsigned j=blocks[i]; j<blocks[i+1]; ++j) { tder.setAtomIndex(n,j); mypos[n]=pos[j]; n++; }
99 453 : 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 : 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 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=std::sqrt(totd); double xx=0.5/totd;
125 15 : myder.scaleAllDerivatives( xx );
126 : }
127 37 : return totd;
128 37 : }
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 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 : 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 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 44 : MultiValue tvals( 1, mypack.getNumberOfDerivatives() ); ReferenceValuePack tder( 0, getNumberOfAtoms(), tvals );
179 132 : for(unsigned i=0; i<domains.size(); ++i) {
180 : // Must extract appropriate positions here
181 88 : 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 88 : if( tder.centeredpos.size()>0 ) tder.rot[0]=mypack.rot[i];
187 : unsigned n=0;
188 396 : for(unsigned j=blocks[i]; j<blocks[i+1]; ++j) {
189 308 : 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 396 : for(unsigned k=n; k<getNumberOfAtoms(); ++k) tder.setAtomIndex(k,3*vecs.size()+10);
197 :
198 : // Do the calculations
199 88 : 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 44 : }
208 :
209 : }
|
