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 : #ifndef __PLUMED_secondarystructure_SecondaryStructureBase_h
23 : #define __PLUMED_secondarystructure_SecondaryStructureBase_h
24 :
25 : #include "core/ActionWithVector.h"
26 : #include "core/ActionShortcut.h"
27 : #include "core/PlumedMain.h"
28 : #include "core/ActionSet.h"
29 : #include "core/GenericMolInfo.h"
30 : #include "core/ParallelTaskManager.h"
31 : #include "tools/ColvarOutput.h"
32 : #include <vector>
33 :
34 : namespace PLMD {
35 : namespace secondarystructure {
36 :
37 : /// Base action for calculating things like AlphRMSD, AntibetaRMSD, etc
38 : template <class T>
39 : class SecondaryStructureBase: public ActionWithVector {
40 : public:
41 : using input_type = T;
42 : using PTM = ParallelTaskManager<SecondaryStructureBase<T>>;
43 : static constexpr size_t virialSize=9;
44 : static constexpr unsigned customGatherStep=3;
45 : static constexpr unsigned customGatherStopBefore=virialSize;
46 :
47 : private:
48 : PTM taskmanager;
49 : public:
50 : static void registerKeywords( Keywords& keys );
51 : static void readBackboneAtoms( ActionShortcut* action, PlumedMain& plumed, const std::string& backnames, std::vector<unsigned>& chain_lengths, std::vector<std::string>& all_atoms );
52 : static bool readShortcutWords( std::string& ltmap, ActionShortcut* action );
53 : explicit SecondaryStructureBase(const ActionOptions&);
54 : unsigned getNumberOfDerivatives() override ;
55 : void calculate() override;
56 : void getInputData( std::vector<double>& inputdata ) const override ;
57 : void applyNonZeroRankForces( std::vector<double>& outforces ) override ;
58 : static void performTask( unsigned task_index, const T& actiondata, ParallelActionsInput& input, ParallelActionsOutput& output );
59 : static int getNumberOfValuesPerTask( std::size_t task_index, const T& actiondata );
60 : static void getForceIndices( std::size_t task_index, std::size_t colno, std::size_t ntotal_force, const T& actiondata, const ParallelActionsInput& input, ForceIndexHolder force_indices );
61 : };
62 :
63 : template <class T>
64 132 : unsigned SecondaryStructureBase<T>::getNumberOfDerivatives() {
65 132 : return 3*getNumberOfAtoms()+virialSize;
66 : }
67 :
68 : template <class T>
69 38 : bool SecondaryStructureBase<T>::readShortcutWords( std::string& ltmap, ActionShortcut* action ) {
70 76 : action->parse("LESS_THAN",ltmap);
71 38 : if( ltmap.length()==0 ) {
72 : std::string nn, mm, d_0, r_0;
73 6 : action->parse("R_0",r_0);
74 3 : if( r_0.length()==0 ) {
75 : r_0="0.08";
76 : }
77 3 : action->parse("NN",nn);
78 3 : action->parse("D_0",d_0);
79 3 : action->parse("MM",mm);
80 6 : ltmap = "RATIONAL R_0=" + r_0 + " D_0=" + d_0 + " NN=" + nn + " MM=" + mm;
81 : return false;
82 : }
83 : return true;
84 : }
85 :
86 : template <class T>
87 288 : void SecondaryStructureBase<T>::registerKeywords( Keywords& keys ) {
88 288 : ActionWithVector::registerKeywords( keys );
89 288 : PTM::registerKeywords( keys );
90 288 : keys.addFlag("NOPBC",false,"ignore the periodic boundary conditions");
91 576 : keys.addInputKeyword("optional","MASK","vector","a vector which is used to determine which elements of the secondary structure variable should be computed");
92 288 : keys.add("atoms","ATOMS","this is the full list of atoms that we are investigating");
93 288 : keys.add("numbered","SEGMENT","this is the lists of atoms in the segment that are being considered");
94 576 : if( keys.getDisplayName()=="SECONDARY_STRUCTURE_DRMSD" ) {
95 36 : keys.add("compulsory","BONDLENGTH","0.17","the length to use for bonds");
96 : }
97 288 : keys.add("numbered","STRUCTURE","the reference structure");
98 576 : if( keys.getDisplayName()=="SECONDARY_STRUCTURE_RMSD" ) {
99 44 : keys.add("compulsory","TYPE","OPTIMAL","the manner in which RMSD alignment is performed. Should be OPTIMAL or SIMPLE. "
100 : "For more details on the OPTIMAL and SIMPLE methods see \\ref RMSD.");
101 488 : } else if( keys.getDisplayName()!="SECONDARY_STRUCTURE_DRMSD" ) {
102 208 : keys.add("compulsory","TYPE","DRMSD","the manner in which RMSD alignment is performed. Should be OPTIMAL, SIMPLE or DRMSD. "
103 : "For more details on the OPTIMAL and SIMPLE methods see \\ref RMSD. For more details on the "
104 : "DRMSD method see \\ref DRMSD.");
105 : }
106 288 : keys.addFlag("VERBOSE",false,"write a more detailed output");
107 576 : keys.reset_style("VERBOSE","hidden");
108 1080 : if( keys.getDisplayName()!="SECONDARY_STRUCTURE_DRMSD" && keys.getDisplayName()!="SECONDARY_STRUCTURE_RMSD" ) {
109 208 : keys.add("residues","RESIDUES","this command is used to specify the set of residues that could conceivably form part of the secondary structure. "
110 : "It is possible to use residues numbers as the various chains and residues should have been identified else using an instance of the "
111 : "\\ref MOLINFO action. If you wish to use all the residues from all the chains in your system you can do so by "
112 : "specifying all. Alternatively, if you wish to use a subset of the residues you can specify the particular residues "
113 : "you are interested in as a list of numbers. Please be aware that to form secondary structure elements your chain "
114 : "must contain at least N residues, where N is dependent on the particular secondary structure you are interested in. "
115 : "As such if you define portions of the chain with fewer than N residues the code will crash.");
116 208 : keys.add("optional","LESS_THAN","calculate the number of a residue segments that are within a certain target distance of this secondary structure type. "
117 : "This quantity is calculated using \\f$\\sum_i \\sigma(s_i)\\f$, where \\f$\\sigma(s)\\f$ is a \\ref switchingfunction.");
118 208 : keys.add("optional","R_0","The r_0 parameter of the switching function.");
119 208 : keys.add("compulsory","D_0","0.0","The d_0 parameter of the switching function");
120 208 : keys.add("compulsory","NN","8","The n parameter of the switching function");
121 208 : keys.add("compulsory","MM","12","The m parameter of the switching function");
122 : }
123 288 : keys.addFlag("ALIGN_STRANDS",false,"ensure that the two halves of a beta sheet are not broken by the periodic boundaries before doing alignment");
124 576 : keys.addOutputComponent("struct","default","scalar","the vectors containing the rmsd distances between the residues and each of the reference structures");
125 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");
126 288 : keys.needsAction("SECONDARY_STRUCTURE_RMSD");
127 288 : keys.needsAction("SECONDARY_STRUCTURE_DRMSD");
128 288 : keys.needsAction("LESS_THAN");
129 288 : keys.needsAction("SUM");
130 288 : keys.addDOI("10.1021/ct900202f");
131 288 : }
132 :
133 : template <class T>
134 38 : void SecondaryStructureBase<T>::readBackboneAtoms( ActionShortcut* action, PlumedMain& plumed, const std::string& moltype, std::vector<unsigned>& chain_lengths, std::vector<std::string>& all_atoms ) {
135 38 : auto* moldat=plumed.getActionSet().selectLatest<GenericMolInfo*>(action);
136 38 : if( ! moldat ) {
137 0 : action->error("Unable to find MOLINFO in input");
138 : }
139 :
140 : std::vector<std::string> resstrings;
141 76 : action->parseVector( "RESIDUES", resstrings );
142 38 : if(resstrings.size()==0) {
143 0 : action->error("residues are not defined, check the keyword RESIDUES");
144 38 : } else if( Tools::caseInSensStringCompare(resstrings[0], "all") ) {
145 : resstrings[0]="all";
146 36 : action->log.printf(" examining all possible secondary structure combinations\n");
147 : } else {
148 2 : action->log.printf(" examining secondary structure in residue positions : %s ",resstrings[0].c_str() );
149 3 : for(unsigned i=1; i<resstrings.size(); ++i) {
150 1 : action->log.printf(", %s",resstrings[i].c_str() );
151 : }
152 2 : action->log.printf("\n");
153 : }
154 : std::vector< std::vector<AtomNumber> > backatoms;
155 38 : moldat->getBackbone( resstrings, moltype, backatoms );
156 :
157 38 : chain_lengths.resize( backatoms.size() );
158 587 : for(unsigned i=0; i<backatoms.size(); ++i) {
159 549 : chain_lengths[i]=backatoms[i].size();
160 22569 : for(unsigned j=0; j<backatoms[i].size(); ++j) {
161 : std::string bat_str;
162 22020 : Tools::convert( backatoms[i][j].serial(), bat_str );
163 22020 : all_atoms.push_back( bat_str );
164 : }
165 : }
166 38 : }
167 :
168 : template <class T>
169 38 : SecondaryStructureBase<T>::SecondaryStructureBase(const ActionOptions&ao):
170 : Action(ao),
171 : ActionWithVector(ao),
172 38 : taskmanager(this) {
173 38 : if( plumed.usingNaturalUnits() ) {
174 0 : error("cannot use this collective variable when using natural units");
175 : }
176 :
177 38 : input_type myinput;
178 38 : parseFlag("NOPBC",myinput.nopbc);
179 38 : std::string alignType="";
180 38 : if( getName()=="SECONDARY_STRUCTURE_RMSD" ) {
181 42 : parse("TYPE",alignType);
182 : }
183 38 : log.printf(" distances from secondary structure elements are calculated using %s algorithm\n", alignType.c_str() );
184 76 : log<<" Bibliography "<<plumed.cite("Pietrucci and Laio, J. Chem. Theory Comput. 5, 2197 (2009)");
185 38 : log<<"\n";
186 :
187 38 : parseFlag("ALIGN_STRANDS",myinput.align_strands);
188 38 : bool verbose_output=false;
189 38 : parseFlag("VERBOSE",verbose_output);
190 38 : log.printf(" ensuring atoms 7 and 22 in each residue are not separated by pbc before doing alignment\n");
191 :
192 : // Read in the atoms
193 : std::vector<AtomNumber> all_atoms;
194 76 : parseAtomList("ATOMS",all_atoms);
195 38 : if( all_atoms.size()==0 ) {
196 0 : error("no atoms were specified -- use ATOMS");
197 : }
198 38 : requestAtoms( all_atoms );
199 :
200 : std::vector<std::vector<unsigned> > colvar_atoms;
201 160063 : for(unsigned i=1;; ++i) {
202 : std::vector<unsigned> newatoms;
203 320202 : if( !parseNumberedVector("SEGMENT",i,newatoms) ) {
204 : break;
205 : }
206 160063 : if( verbose_output ) {
207 0 : log.printf(" Secondary structure segment %u contains atoms : ", static_cast<unsigned>(colvar_atoms.size()+1));
208 0 : for(unsigned ii=0; ii<newatoms.size(); ++ii) {
209 0 : log.printf("%d ",all_atoms[newatoms[ii]].serial() );
210 : }
211 0 : log.printf("\n");
212 : }
213 160063 : colvar_atoms.push_back( newatoms );
214 : }
215 38 : if( colvar_atoms.size()==0 ) {
216 0 : error("did not find any SEGMENT keywords in input");
217 : }
218 38 : myinput.colvar_atoms.resize( colvar_atoms.size(), colvar_atoms[0].size() );
219 160101 : for(unsigned i=0; i<colvar_atoms.size(); ++i) {
220 4961953 : for(unsigned j=0; j<colvar_atoms[i].size(); ++j) {
221 4801890 : myinput.colvar_atoms[i][j] = colvar_atoms[i][j];
222 : }
223 : }
224 :
225 38 : double bondlength=0.0;
226 38 : if( getName()=="SECONDARY_STRUCTURE_DRMSD" ) {
227 17 : parse("BONDLENGTH",bondlength);
228 17 : bondlength=bondlength/getUnits().getLength();
229 : }
230 : // Read in the reference structure
231 56 : for(unsigned ii=1;; ++ii) {
232 : std::vector<double> cstruct;
233 188 : if( !parseNumberedVector("STRUCTURE",ii,cstruct) ) {
234 : break ;
235 : }
236 56 : plumed_assert( cstruct.size()%3==0 && cstruct.size()/3==colvar_atoms[0].size() );
237 56 : std::vector<Vector> structure( cstruct.size()/3 );
238 1736 : for(unsigned i=0; i<structure.size(); ++i) {
239 6720 : for(unsigned j=0; j<3; ++j) {
240 5040 : structure[i][j] = 0.1*cstruct[3*i+j]/getUnits().getLength();
241 : }
242 : }
243 81 : myinput.setReferenceStructure( alignType, bondlength, structure );
244 : }
245 :
246 : // And create values to hold everything
247 38 : plumed_assert( myinput.nstructures>0 );
248 38 : std::vector<std::size_t> shape(1);
249 38 : shape[0]=colvar_atoms.size();
250 38 : if( myinput.nstructures==1 ) {
251 20 : addValue( shape );
252 20 : setNotPeriodic();
253 : } else {
254 : std::string num;
255 54 : for(unsigned i=0; i<myinput.nstructures; ++i) {
256 36 : Tools::convert( i+1, num );
257 36 : addComponent( "struct-" + num, shape );
258 72 : componentIsNotPeriodic( "struct-" + num );
259 : }
260 : }
261 94 : for(unsigned i=0; i<getNumberOfComponents(); ++i) {
262 56 : getPntrToComponent(i)->setDerivativeIsZeroWhenValueIsZero();
263 : }
264 38 : myinput.nindices_per_task = colvar_atoms[0].size();
265 38 : taskmanager.setupParallelTaskManager( 3*colvar_atoms[0].size() + virialSize, 3*getNumberOfAtoms() + 9 );
266 38 : taskmanager.setActionInput( myinput );
267 76 : }
268 :
269 : template <class T>
270 276 : void SecondaryStructureBase<T>::calculate() {
271 276 : taskmanager.runAllTasks();
272 276 : }
273 :
274 : template <class T>
275 276 : void SecondaryStructureBase<T>::getInputData( std::vector<double>& inputdata ) const {
276 276 : if( inputdata.size()!=3*getNumberOfAtoms() ) {
277 38 : inputdata.resize( 3*getNumberOfAtoms() );
278 : }
279 :
280 : std::size_t k=0;
281 134916 : for(unsigned i=0; i<getNumberOfAtoms(); ++i) {
282 134640 : Vector mypos( getPosition(i) );
283 134640 : inputdata[k] = mypos[0];
284 134640 : k++;
285 134640 : inputdata[k] = mypos[1];
286 134640 : k++;
287 134640 : inputdata[k] = mypos[2];
288 134640 : k++;
289 : }
290 276 : }
291 :
292 : template <class T>
293 140516 : void SecondaryStructureBase<T>::performTask( unsigned task_index, const T& actiondata, ParallelActionsInput& input, ParallelActionsOutput& output ) {
294 : // std::vector<Vector> pos( actiondata.natoms );
295 : std::array<Vector,30> pos;
296 :
297 4355996 : for(unsigned i=0; i<actiondata.natoms; ++i) {
298 4215480 : const unsigned atno = actiondata.colvar_atoms(task_index,i);
299 4215480 : pos[i][0] = input.inputdata[3*atno+0];
300 4215480 : pos[i][1] = input.inputdata[3*atno+1];
301 4215480 : pos[i][2] = input.inputdata[3*atno+2];
302 : }
303 : // This aligns the two strands if this is required
304 140516 : if( actiondata.align_strands ) {
305 140240 : Vector distance=input.pbc->distance( pos[6],pos[21] );
306 : Vector origin_old, origin_new;
307 140240 : origin_old=pos[21];
308 : origin_new=pos[6]+distance;
309 2243840 : for(unsigned i=15; i<30; ++i) {
310 2103600 : pos[i]+=( origin_new - origin_old );
311 : }
312 276 : } else if( !actiondata.nopbc ) {
313 8280 : for(unsigned i=0; i<actiondata.natoms-1; ++i) {
314 : const Vector & first (pos[i]);
315 8004 : Vector & second (pos[i+1]);
316 8004 : second=first+input.pbc->distance(first,second);
317 : }
318 : }
319 :
320 : // Create a holder for the derivatives
321 140516 : const unsigned rs = actiondata.nstructures;
322 : ColvarOutput rmsd_output( output.values,
323 : 3*pos.size()+virialSize,
324 : output.derivatives.data() );
325 : // And now calculate the DRMSD
326 354150 : for(unsigned i=0; i<rs; ++i) {
327 213634 : T::calculateDistance( i, input.noderiv, actiondata, View{pos.data(),pos.size()}, rmsd_output );
328 : }
329 140516 : }
330 :
331 : template <class T>
332 240 : void SecondaryStructureBase<T>::applyNonZeroRankForces( std::vector<double>& outforces ) {
333 240 : taskmanager.applyForces( outforces );
334 240 : }
335 :
336 : template <class T>
337 70096 : int SecondaryStructureBase<T>::getNumberOfValuesPerTask( std::size_t task_index,
338 : const T& actiondata ) {
339 70096 : return 1;
340 : }
341 :
342 : template <class T>
343 70096 : void SecondaryStructureBase<T>::getForceIndices( std::size_t task_index,
344 : std::size_t /* colno */,
345 : std::size_t ntotal_force,
346 : const T& actiondata,
347 : const ParallelActionsInput& input,
348 : ForceIndexHolder force_indices ) {
349 176631 : for(unsigned i=0; i<input.ncomponents; ++i) {
350 : std::size_t m = 0;
351 3302585 : for(unsigned j=0; j<actiondata.nindices_per_task; ++j) {
352 3196050 : std::size_t base = 3*actiondata.colvar_atoms[task_index][j];
353 3196050 : force_indices.indices[i][m] = base + 0;
354 3196050 : ++m;
355 3196050 : force_indices.indices[i][m] = base + 1;
356 3196050 : ++m;
357 3196050 : force_indices.indices[i][m] = base + 2;
358 3196050 : ++m;
359 : }
360 1065350 : for(unsigned n=ntotal_force-virialSize; n<ntotal_force; ++n) {
361 958815 : force_indices.indices[i][m] = n;
362 958815 : ++m;
363 : }
364 106535 : force_indices.threadsafe_derivatives_end[i] = 0;
365 106535 : force_indices.tot_indices[i] = m;
366 : }
367 70096 : }
368 :
369 : }
370 : }
371 :
372 : #endif
|
