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