Line data Source code
1 : /* +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
2 : Copyright (c) 2013-2020 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 "SecondaryStructureRMSD.h"
23 : #include "core/PlumedMain.h"
24 : #include "core/ActionSet.h"
25 : #include "core/GenericMolInfo.h"
26 : #include "core/ActionShortcut.h"
27 : #include "core/ActionRegister.h"
28 :
29 : //+PLUMEDOC MCOLVAR SECONDARY_STRUCTURE_RMSD
30 : /*
31 : Calclulate the distance between segments of a protein and a reference structure of interest
32 :
33 : \par Examples
34 :
35 : */
36 : //+ENDPLUMEDOC
37 :
38 : namespace PLMD {
39 : namespace secondarystructure {
40 :
41 : PLUMED_REGISTER_ACTION(SecondaryStructureRMSD,"SECONDARY_STRUCTURE_RMSD");
42 :
43 38 : bool SecondaryStructureRMSD::readShortcutWords( std::string& ltmap, ActionShortcut* action ) {
44 76 : action->parse("LESS_THAN",ltmap);
45 38 : if( ltmap.length()==0 ) {
46 : std::string nn, mm, d_0, r_0;
47 6 : action->parse("R_0",r_0);
48 3 : if( r_0.length()==0 ) {
49 : r_0="0.08";
50 : }
51 3 : action->parse("NN",nn);
52 3 : action->parse("D_0",d_0);
53 3 : action->parse("MM",mm);
54 6 : ltmap = "RATIONAL R_0=" + r_0 + " D_0=" + d_0 + " NN=" + nn + " MM=" + mm;
55 : return false;
56 : }
57 : return true;
58 : }
59 :
60 38 : void SecondaryStructureRMSD::expandShortcut( const bool& uselessthan, const std::string& labout, const std::string& labin, const std::string& ltmap, ActionShortcut* action ) {
61 76 : action->readInputLine( labout + "_lt: LESS_THAN ARG=" + labin + " SWITCH={" + ltmap +"}");
62 38 : if( uselessthan ) {
63 70 : action->readInputLine( labout + "_lessthan: SUM ARG=" + labout + "_lt PERIODIC=NO");
64 : } else {
65 6 : action->readInputLine( labout + ": SUM ARG=" + labout + "_lt PERIODIC=NO");
66 : }
67 38 : }
68 :
69 288 : void SecondaryStructureRMSD::registerKeywords( Keywords& keys ) {
70 288 : ActionWithVector::registerKeywords( keys );
71 576 : keys.addFlag("NOPBC",false,"ignore the periodic boundary conditions");
72 576 : keys.add("residues","RESIDUES","this command is used to specify the set of residues that could conceivably form part of the secondary structure. "
73 : "It is possible to use residues numbers as the various chains and residues should have been identified else using an instance of the "
74 : "\\ref MOLINFO action. If you wish to use all the residues from all the chains in your system you can do so by "
75 : "specifying all. Alternatively, if you wish to use a subset of the residues you can specify the particular residues "
76 : "you are interested in as a list of numbers. Please be aware that to form secondary structure elements your chain "
77 : "must contain at least N residues, where N is dependent on the particular secondary structure you are interested in. "
78 : "As such if you define portions of the chain with fewer than N residues the code will crash.");
79 576 : keys.add("atoms","ATOMS","this is the full list of atoms that we are investigating");
80 576 : keys.add("numbered","SEGMENT","this is the lists of atoms in the segment that are being considered");
81 576 : keys.add("compulsory","BONDLENGTH","the length to use for bonds");
82 576 : keys.add("numbered","STRUCTURE","the reference structure");
83 576 : keys.add("compulsory","TYPE","DRMSD","the manner in which RMSD alignment is performed. Should be OPTIMAL, SIMPLE or DRMSD. "
84 : "For more details on the OPTIMAL and SIMPLE methods see \\ref RMSD. For more details on the "
85 : "DRMSD method see \\ref DRMSD.");
86 576 : keys.add("optional","STRANDS_CUTOFF","If in a segment of protein the two strands are further apart then the calculation "
87 : "of the actual RMSD is skipped as the structure is very far from being beta-sheet like. "
88 : "This keyword speeds up the calculation enormously when you are using the LESS_THAN option. "
89 : "However, if you are using some other option, then this cannot be used");
90 576 : keys.add("optional","CUTOFF_ATOMS","the pair of atoms that are used to calculate the strand cutoff");
91 576 : keys.addFlag("VERBOSE",false,"write a more detailed output");
92 576 : keys.add("optional","LESS_THAN","calculate the number of a residue segments that are within a certain target distance of this secondary structure type. "
93 : "This quantity is calculated using \\f$\\sum_i \\sigma(s_i)\\f$, where \\f$\\sigma(s)\\f$ is a \\ref switchingfunction.");
94 576 : keys.add("optional","R_0","The r_0 parameter of the switching function.");
95 576 : keys.add("compulsory","D_0","0.0","The d_0 parameter of the switching function");
96 576 : keys.add("compulsory","NN","8","The n parameter of the switching function");
97 576 : keys.add("compulsory","MM","12","The m parameter of the switching function");
98 576 : keys.add("hidden","NO_ACTION_LOG","suppresses printing from action on the log");
99 576 : keys.addOutputComponent("struct","default","the vectors containing the rmsd distances between the residues and each of the reference structures");
100 576 : keys.addOutputComponent("lessthan","default","the number blocks of residues that have an RMSD from the secondary structure that is less than the threshold");
101 288 : keys.needsAction("SECONDARY_STRUCTURE_RMSD");
102 288 : keys.needsAction("LESS_THAN");
103 288 : keys.needsAction("SUM");
104 288 : }
105 :
106 38 : void SecondaryStructureRMSD::readBackboneAtoms( ActionShortcut* action, PlumedMain& plumed, const std::string& moltype, std::vector<unsigned>& chain_lengths, std::string& all_atoms ) {
107 38 : auto* moldat=plumed.getActionSet().selectLatest<GenericMolInfo*>(action);
108 38 : if( ! moldat ) {
109 0 : action->error("Unable to find MOLINFO in input");
110 : }
111 :
112 : std::vector<std::string> resstrings;
113 76 : action->parseVector( "RESIDUES", resstrings );
114 38 : if(resstrings.size()==0) {
115 0 : action->error("residues are not defined, check the keyword RESIDUES");
116 76 : } else if( Tools::caseInSensStringCompare(resstrings[0], "all") ) {
117 : resstrings[0]="all";
118 36 : action->log.printf(" examining all possible secondary structure combinations\n");
119 : } else {
120 2 : action->log.printf(" examining secondary structure in residue positions : %s ",resstrings[0].c_str() );
121 3 : for(unsigned i=1; i<resstrings.size(); ++i) {
122 1 : action->log.printf(", %s",resstrings[i].c_str() );
123 : }
124 2 : action->log.printf("\n");
125 : }
126 : std::vector< std::vector<AtomNumber> > backatoms;
127 38 : moldat->getBackbone( resstrings, moltype, backatoms );
128 :
129 38 : chain_lengths.resize( backatoms.size() );
130 587 : for(unsigned i=0; i<backatoms.size(); ++i) {
131 549 : chain_lengths[i]=backatoms[i].size();
132 22569 : for(unsigned j=0; j<backatoms[i].size(); ++j) {
133 : std::string bat_str;
134 22020 : Tools::convert( backatoms[i][j].serial(), bat_str );
135 22020 : if( i==0 && j==0 ) {
136 76 : all_atoms = "ATOMS=" + bat_str;
137 : } else {
138 43964 : all_atoms += "," + bat_str;
139 : }
140 : }
141 : }
142 38 : }
143 :
144 :
145 38 : SecondaryStructureRMSD::SecondaryStructureRMSD(const ActionOptions&ao):
146 : Action(ao),
147 : ActionWithVector(ao),
148 38 : nopbc(false),
149 38 : align_strands(false),
150 38 : s_cutoff2(0),
151 38 : align_atom_1(0),
152 38 : align_atom_2(0) {
153 38 : if( plumed.usingNaturalUnits() ) {
154 0 : error("cannot use this collective variable when using natural units");
155 : }
156 :
157 38 : parse("TYPE",alignType);
158 38 : parseFlag("NOPBC",nopbc);
159 38 : log.printf(" distances from secondary structure elements are calculated using %s algorithm\n",alignType.c_str() );
160 76 : log<<" Bibliography "<<plumed.cite("Pietrucci and Laio, J. Chem. Theory Comput. 5, 2197 (2009)");
161 38 : log<<"\n";
162 :
163 38 : parseFlag("VERBOSE",verbose_output);
164 :
165 76 : if( keywords.exists("STRANDS_CUTOFF") ) {
166 38 : double s_cutoff = 0;
167 38 : parse("STRANDS_CUTOFF",s_cutoff);
168 38 : align_strands=true;
169 38 : if( s_cutoff>0) {
170 32 : log.printf(" ignoring contributions from strands that are more than %f apart\n",s_cutoff);
171 : std::vector<unsigned> cutatoms;
172 64 : parseVector("CUTOFF_ATOMS",cutatoms);
173 32 : if( cutatoms.size()==2 ) {
174 32 : align_atom_1=cutatoms[0];
175 32 : align_atom_2=cutatoms[1];
176 : } else {
177 0 : error("did not find CUTOFF_ATOMS in input");
178 : }
179 : }
180 38 : s_cutoff2=s_cutoff*s_cutoff;
181 : }
182 :
183 : // Read in the atoms
184 : std::vector<AtomNumber> all_atoms;
185 38 : parseAtomList("ATOMS",all_atoms);
186 38 : requestAtoms( all_atoms );
187 :
188 160063 : for(unsigned i=1;; ++i) {
189 : std::vector<unsigned> newatoms;
190 320202 : if( !parseNumberedVector("SEGMENT",i,newatoms) ) {
191 : break;
192 : }
193 160063 : if( verbose_output ) {
194 0 : log.printf(" Secondary structure segment %u contains atoms : ", static_cast<unsigned>(colvar_atoms.size()+1));
195 0 : for(unsigned i=0; i<newatoms.size(); ++i) {
196 0 : log.printf("%d ",all_atoms[newatoms[i]].serial() );
197 : }
198 0 : log.printf("\n");
199 : }
200 160063 : colvar_atoms.push_back( newatoms );
201 160063 : }
202 38 : if( colvar_atoms.size()==0 ) {
203 0 : error("missing SEGMENT keywords -- no atoms have been specified for comparison");
204 : }
205 :
206 : double bondlength;
207 38 : parse("BONDLENGTH",bondlength);
208 38 : bondlength=bondlength/getUnits().getLength();
209 :
210 : // Read in the reference structure
211 56 : for(unsigned ii=1;; ++ii) {
212 : std::vector<double> cstruct;
213 188 : if( !parseNumberedVector("STRUCTURE",ii,cstruct) ) {
214 : break ;
215 : }
216 56 : plumed_assert( cstruct.size()%3==0 && cstruct.size()/3==colvar_atoms[0].size() );
217 56 : std::vector<Vector> structure( cstruct.size()/3 );
218 1736 : for(unsigned i=0; i<structure.size(); ++i) {
219 6720 : for(unsigned j=0; j<3; ++j) {
220 5040 : structure[i][j] = 0.1*cstruct[3*i+j]/getUnits().getLength();
221 : }
222 : }
223 56 : if( alignType=="DRMSD" ) {
224 : std::map<std::pair<unsigned,unsigned>, double> targets;
225 750 : for(unsigned i=0; i<structure.size()-1; ++i) {
226 11600 : for(unsigned j=i+1; j<structure.size(); ++j) {
227 10875 : double distance = delta( structure[i], structure[j] ).modulo();
228 10875 : if(distance > bondlength) {
229 10172 : targets[std::make_pair(i,j)] = distance;
230 : }
231 : }
232 : }
233 25 : drmsd_targets.push_back( targets );
234 : } else {
235 31 : Vector center;
236 31 : std::vector<double> align( structure.size(), 1.0 ), displace( structure.size(), 1.0 );
237 961 : for(unsigned i=0; i<structure.size(); ++i) {
238 930 : center+=structure[i]*align[i];
239 : }
240 961 : for(unsigned i=0; i<structure.size(); ++i) {
241 930 : structure[i] -= center;
242 : }
243 31 : RMSD newrmsd;
244 31 : newrmsd.clear();
245 31 : newrmsd.set(align,displace,structure,alignType,true,true);
246 31 : myrmsd.push_back( newrmsd );
247 31 : }
248 56 : }
249 :
250 : // And create values to hold everything
251 : unsigned nref = myrmsd.size();
252 38 : if( alignType=="DRMSD" ) {
253 : nref=drmsd_targets.size();
254 : }
255 38 : plumed_assert( nref>0 );
256 38 : std::vector<unsigned> shape(1);
257 38 : shape[0]=colvar_atoms.size();
258 38 : if( nref==1 ) {
259 20 : addValue( shape );
260 20 : setNotPeriodic();
261 : } else {
262 : std::string num;
263 54 : for(unsigned i=0; i<nref; ++i) {
264 36 : Tools::convert( i+1, num );
265 36 : addComponent( "struct-" + num, shape );
266 72 : componentIsNotPeriodic( "struct-" + num );
267 : }
268 : }
269 38 : }
270 :
271 94 : void SecondaryStructureRMSD::areAllTasksRequired( std::vector<ActionWithVector*>& task_reducing_actions ) {
272 94 : if( s_cutoff2>0 ) {
273 80 : task_reducing_actions.push_back(this);
274 : }
275 94 : }
276 :
277 761700 : int SecondaryStructureRMSD::checkTaskStatus( const unsigned& taskno, int& flag ) const {
278 761700 : if( s_cutoff2>0 ) {
279 761700 : Vector distance=pbcDistance( ActionAtomistic::getPosition( getAtomIndex(taskno,align_atom_1) ),
280 : ActionAtomistic::getPosition( getAtomIndex(taskno,align_atom_2) ) );
281 761700 : if( distance.modulo2()<s_cutoff2 ) {
282 : return 1;
283 : }
284 691712 : return 0;
285 : }
286 0 : return flag;
287 : }
288 :
289 276 : void SecondaryStructureRMSD::calculate() {
290 276 : runAllTasks();
291 276 : }
292 :
293 70420 : void SecondaryStructureRMSD::performTask( const unsigned& current, MultiValue& myvals ) const {
294 : // Resize the derivatives if need be
295 70420 : unsigned nderi = 3*getNumberOfAtoms()+9;
296 70420 : if( myvals.getNumberOfDerivatives()!=nderi ) {
297 36 : myvals.resize( myvals.getNumberOfValues(), nderi, 0, 0 );
298 : }
299 : // Retrieve the positions
300 70420 : const unsigned natoms = colvar_atoms[current].size();
301 70420 : std::vector<Vector> pos( natoms ), deriv( natoms );
302 2183020 : for(unsigned i=0; i<natoms; ++i) {
303 2112600 : pos[i]=ActionAtomistic::getPosition( getAtomIndex(current,i) );
304 : }
305 :
306 : // This aligns the two strands if this is required
307 70420 : Vector distance=pbcDistance( pos[align_atom_1],pos[align_atom_2] );
308 70420 : if( align_strands ) {
309 70420 : Vector origin_old, origin_new;
310 70420 : origin_old=pos[align_atom_2];
311 70420 : origin_new=pos[align_atom_1]+distance;
312 1126720 : for(unsigned i=15; i<30; ++i) {
313 1056300 : pos[i]+=( origin_new - origin_old );
314 : }
315 0 : } else if( !nopbc ) {
316 0 : for(unsigned i=0; i<natoms-1; ++i) {
317 0 : const Vector & first (pos[i]);
318 0 : Vector & second (pos[i+1]);
319 0 : second=first+pbcDistance(first,second);
320 : }
321 : }
322 : // Create a holder for the derivatives
323 70420 : if( alignType=="DRMSD" ) {
324 : // And now calculate the DRMSD
325 21864 : const unsigned rs = drmsd_targets.size();
326 56167 : for(unsigned i=0; i<rs; ++i) {
327 : double drmsd=0;
328 34303 : Vector distance;
329 34303 : Tensor vir;
330 34303 : vir.zero();
331 1063393 : for(unsigned j=0; j<natoms; ++j) {
332 1029090 : deriv[j].zero();
333 : }
334 13995444 : for(const auto & it : drmsd_targets[i] ) {
335 13961141 : const unsigned k=it.first.first;
336 13961141 : const unsigned j=it.first.second;
337 :
338 13961141 : distance=delta( pos[k], pos[j] );
339 13961141 : const double len = distance.modulo();
340 13961141 : const double diff = len - it.second;
341 13961141 : const double der = diff / len;
342 13961141 : drmsd += diff*diff;
343 :
344 13961141 : if( !doNotCalculateDerivatives() ) {
345 13731677 : deriv[k] += -der*distance;
346 13731677 : deriv[j] += der*distance;
347 13731677 : vir += -der*Tensor(distance,distance);
348 : }
349 : }
350 :
351 34303 : const double inpairs = 1./static_cast<double>(drmsd_targets[i].size());
352 34303 : unsigned ostrn = getConstPntrToComponent(i)->getPositionInStream();
353 34303 : drmsd = sqrt(inpairs*drmsd);
354 34303 : myvals.setValue( ostrn, drmsd );
355 :
356 34303 : if( !doNotCalculateDerivatives() ) {
357 33739 : double scalef = inpairs / drmsd;
358 1045909 : for(unsigned j=0; j<natoms; ++j) {
359 : const unsigned ja = getAtomIndex( current, j );
360 1012170 : myvals.addDerivative( ostrn, 3*ja + 0, scalef*deriv[j][0] );
361 1012170 : myvals.updateIndex( ostrn, 3*ja+0 );
362 1012170 : myvals.addDerivative( ostrn, 3*ja + 1, scalef*deriv[j][1] );
363 1012170 : myvals.updateIndex( ostrn, 3*ja+1 );
364 1012170 : myvals.addDerivative( ostrn, 3*ja + 2, scalef*deriv[j][2] );
365 1012170 : myvals.updateIndex( ostrn, 3*ja+2 );
366 : }
367 33739 : unsigned nbase = myvals.getNumberOfDerivatives() - 9;
368 134956 : for(unsigned k=0; k<3; ++k) {
369 404868 : for(unsigned j=0; j<3; ++j) {
370 303651 : myvals.addDerivative( ostrn, nbase + 3*k + j, scalef*vir(k,j) );
371 303651 : myvals.updateIndex( ostrn, nbase + 3*k + j );
372 : }
373 : }
374 : }
375 : }
376 : } else {
377 48556 : const unsigned rs = myrmsd.size();
378 121352 : for(unsigned i=0; i<rs; ++i) {
379 72796 : double nr = myrmsd[i].calculate( pos, deriv, false );
380 72796 : unsigned ostrn = getConstPntrToComponent(i)->getPositionInStream();
381 72796 : myvals.setValue( ostrn, nr );
382 :
383 72796 : if( !doNotCalculateDerivatives() ) {
384 72796 : Tensor vir;
385 72796 : vir.zero();
386 2256676 : for(unsigned j=0; j<natoms; ++j) {
387 : const unsigned ja = getAtomIndex( current, j );
388 2183880 : myvals.addDerivative( ostrn, 3*ja + 0, deriv[j][0] );
389 2183880 : myvals.updateIndex( ostrn, 3*colvar_atoms[current][j]+0 );
390 2183880 : myvals.addDerivative( ostrn, 3*ja + 1, deriv[j][1] );
391 2183880 : myvals.updateIndex( ostrn, 3*colvar_atoms[current][j]+1 );
392 2183880 : myvals.addDerivative( ostrn, 3*ja + 2, deriv[j][2] );
393 2183880 : myvals.updateIndex( ostrn, 3*colvar_atoms[current][j]+2 );
394 2183880 : vir+=(-1.0*Tensor( pos[j], deriv[j] ));
395 : }
396 72796 : unsigned nbase = myvals.getNumberOfDerivatives() - 9;
397 291184 : for(unsigned k=0; k<3; ++k) {
398 873552 : for(unsigned j=0; j<3; ++j) {
399 655164 : myvals.addDerivative( ostrn, nbase + 3*k + j, vir(k,j) );
400 655164 : myvals.updateIndex( ostrn, nbase + 3*k + j );
401 : }
402 : }
403 : }
404 : }
405 : }
406 70420 : return;
407 : }
408 :
409 : }
410 : }
|
