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 "SecondaryStructureRMSD.h"
23 : #include "core/PlumedMain.h"
24 : #include "core/ActionSet.h"
25 : #include "core/GenericMolInfo.h"
26 : #include "core/Atoms.h"
27 : #include "vesselbase/Vessel.h"
28 : #include "reference/MetricRegister.h"
29 : #include "reference/SingleDomainRMSD.h"
30 :
31 : namespace PLMD {
32 : namespace secondarystructure {
33 :
34 55 : void SecondaryStructureRMSD::registerKeywords( Keywords& keys ) {
35 55 : Action::registerKeywords( keys );
36 55 : ActionWithValue::registerKeywords( keys );
37 55 : ActionAtomistic::registerKeywords( keys );
38 110 : keys.add("residues","RESIDUES","this command is used to specify the set of residues that could conceivably form part of the secondary structure. "
39 : "It is possible to use residues numbers as the various chains and residues should have been identified else using an instance of the "
40 : "\\ref MOLINFO action. If you wish to use all the residues from all the chains in your system you can do so by "
41 : "specifying all. Alternatively, if you wish to use a subset of the residues you can specify the particular residues "
42 : "you are interested in as a list of numbers. Please be aware that to form secondary structure elements your chain "
43 : "must contain at least N residues, where N is dependent on the particular secondary structure you are interested in. "
44 : "As such if you define portions of the chain with fewer than N residues the code will crash.");
45 110 : keys.add("compulsory","TYPE","DRMSD","the manner in which RMSD alignment is performed. Should be OPTIMAL, SIMPLE or DRMSD. "
46 : "For more details on the OPTIMAL and SIMPLE methods see \\ref RMSD. For more details on the "
47 : "DRMSD method see \\ref DRMSD.");
48 110 : keys.addFlag("NOPBC",false,"ignore the periodic boundary conditions");
49 110 : keys.add("compulsory","R_0","0.08","The r_0 parameter of the switching function.");
50 110 : keys.add("compulsory","D_0","0.0","The d_0 parameter of the switching function");
51 110 : keys.add("compulsory","NN","8","The n parameter of the switching function");
52 110 : keys.add("compulsory","MM","12","The m parameter of the switching function");
53 110 : keys.reserve("optional","STRANDS_CUTOFF","If in a segment of protein the two strands are further apart then the calculation "
54 : "of the actual RMSD is skipped as the structure is very far from being beta-sheet like. "
55 : "This keyword speeds up the calculation enormously when you are using the LESS_THAN option. "
56 : "However, if you are using some other option, then this cannot be used");
57 110 : keys.addFlag("VERBOSE",false,"write a more detailed output");
58 110 : keys.add("hidden","NL_STRIDE","the frequency with which the neighbor list should be updated. Between neighbour list update steps all quantities "
59 : "that contributed less than TOL at the previous neighbor list update step are ignored.");
60 55 : ActionWithVessel::registerKeywords( keys );
61 55 : keys.use("LESS_THAN");
62 55 : keys.use("MIN");
63 55 : keys.use("ALT_MIN");
64 55 : keys.use("LOWEST");
65 55 : keys.use("HIGHEST");
66 55 : keys.setComponentsIntroduction("By default this Action calculates the number of structural units that are within a certain "
67 : "distance of a idealized secondary structure element. This quantity can then be referenced "
68 : "elsewhere in the input by using the label of the action. However, this Action can also be used to "
69 : "calculate the following quantities by using the keywords as described below. The quantities then "
70 : "calculated can be referenced using the label of the action followed by a dot and then the name "
71 : "from the table below. Please note that you can use the LESS_THAN keyword more than once. The resulting "
72 : "components will be labelled <em>label</em>.lessthan-1, <em>label</em>.lessthan-2 and so on unless you "
73 : "exploit the fact that these labels can be given custom labels by using the LABEL keyword in the "
74 : "description of you LESS_THAN function that you are computing");
75 55 : }
76 :
77 43 : SecondaryStructureRMSD::SecondaryStructureRMSD(const ActionOptions&ao):
78 : Action(ao),
79 : ActionAtomistic(ao),
80 : ActionWithValue(ao),
81 : ActionWithVessel(ao),
82 43 : nopbc(false),
83 43 : align_strands(false),
84 43 : s_cutoff2(0),
85 43 : align_atom_1(0),
86 43 : align_atom_2(0) {
87 43 : parse("TYPE",alignType);
88 43 : parseFlag("NOPBC",nopbc);
89 43 : log.printf(" distances from secondary structure elements are calculated using %s algorithm\n",alignType.c_str() );
90 86 : log<<" Bibliography "<<plumed.cite("Pietrucci and Laio, J. Chem. Theory Comput. 5, 2197 (2009)");
91 43 : log<<"\n";
92 :
93 43 : parseFlag("VERBOSE",verbose_output);
94 :
95 86 : if( keywords.exists("STRANDS_CUTOFF") ) {
96 40 : double s_cutoff = 0;
97 40 : parse("STRANDS_CUTOFF",s_cutoff);
98 40 : align_strands=true;
99 40 : if( s_cutoff>0) {
100 38 : log.printf(" ignoring contributions from strands that are more than %f apart\n",s_cutoff);
101 : }
102 40 : s_cutoff2=s_cutoff*s_cutoff;
103 : }
104 43 : }
105 :
106 43 : SecondaryStructureRMSD::~SecondaryStructureRMSD() {
107 : // destructor needed to delete forward declarated objects
108 86 : }
109 :
110 110 : void SecondaryStructureRMSD::turnOnDerivatives() {
111 110 : ActionWithValue::turnOnDerivatives();
112 110 : needsDerivatives();
113 110 : }
114 :
115 40 : void SecondaryStructureRMSD::setAtomsFromStrands( const unsigned& atom1, const unsigned& atom2 ) {
116 40 : align_atom_1=atom1;
117 40 : align_atom_2=atom2;
118 40 : }
119 :
120 43 : void SecondaryStructureRMSD::readBackboneAtoms( const std::string& moltype, std::vector<unsigned>& chain_lengths ) {
121 43 : auto* moldat=plumed.getActionSet().selectLatest<GenericMolInfo*>(this);
122 43 : if( ! moldat ) {
123 0 : error("Unable to find MOLINFO in input");
124 : }
125 :
126 : std::vector<std::string> resstrings;
127 43 : parseVector( "RESIDUES", resstrings );
128 43 : if( !verbose_output ) {
129 43 : if(resstrings.size()==0) {
130 0 : error("residues are not defined, check the keyword RESIDUES");
131 43 : } else if(resstrings[0]=="all") {
132 39 : log.printf(" examining all possible secondary structure combinations\n");
133 : } else {
134 4 : log.printf(" examining secondary structure in residue positions : %s \n",resstrings[0].c_str() );
135 6 : for(unsigned i=1; i<resstrings.size(); ++i) {
136 2 : log.printf(", %s",resstrings[i].c_str() );
137 : }
138 4 : log.printf("\n");
139 : }
140 : }
141 : std::vector< std::vector<AtomNumber> > backatoms;
142 43 : moldat->getBackbone( resstrings, moltype, backatoms );
143 :
144 43 : chain_lengths.resize( backatoms.size() );
145 700 : for(unsigned i=0; i<backatoms.size(); ++i) {
146 657 : chain_lengths[i]=backatoms[i].size();
147 26877 : for(unsigned j=0; j<backatoms[i].size(); ++j) {
148 26220 : all_atoms.push_back( backatoms[i][j] );
149 : }
150 : }
151 43 : ActionAtomistic::requestAtoms( all_atoms );
152 43 : forcesToApply.resize( getNumberOfDerivatives() );
153 43 : }
154 :
155 187632 : void SecondaryStructureRMSD::addColvar( const std::vector<unsigned>& newatoms ) {
156 187632 : if( colvar_atoms.size()>0 ) {
157 187590 : plumed_assert( colvar_atoms[0].size()==newatoms.size() );
158 : }
159 187632 : if( verbose_output ) {
160 0 : log.printf(" Secondary structure segment %u contains atoms : ", static_cast<unsigned>(colvar_atoms.size()+1));
161 0 : for(unsigned i=0; i<newatoms.size(); ++i) {
162 0 : log.printf("%d ",all_atoms[newatoms[i]].serial() );
163 : }
164 0 : log.printf("\n");
165 : }
166 187632 : addTaskToList( colvar_atoms.size() );
167 187632 : colvar_atoms.push_back( newatoms );
168 187632 : }
169 :
170 62 : void SecondaryStructureRMSD::setSecondaryStructure( std::vector<Vector>& structure, double bondlength, double units ) {
171 : // If we are in natural units get conversion factor from nm into natural length units
172 62 : if( plumed.getAtoms().usingNaturalUnits() ) {
173 0 : error("cannot use this collective variable when using natural units");
174 : }
175 62 : plumed_massert( !(align_strands && align_atom_1==0 && align_atom_2==0), "you must use setAtomsFromStrands with strands cutoff");
176 :
177 : // Convert into correct units
178 1922 : for(unsigned i=0; i<structure.size(); ++i) {
179 1860 : structure[i][0]*=units;
180 1860 : structure[i][1]*=units;
181 1860 : structure[i][2]*=units;
182 : }
183 :
184 62 : if( references.size()==0 ) {
185 43 : readVesselKeywords();
186 43 : if( getNumberOfVessels()==0 ) {
187 : double r0;
188 2 : parse("R_0",r0);
189 : double d0;
190 2 : parse("D_0",d0);
191 : int nn;
192 2 : parse("NN",nn);
193 : int mm;
194 2 : parse("MM",mm);
195 2 : std::ostringstream ostr;
196 2 : ostr<<"RATIONAL R_0="<<r0<<" D_0="<<d0<<" NN="<<nn<<" MM="<<mm;
197 : std::string input=ostr.str();
198 2 : addVessel( "LESS_THAN", input, -1 ); // -1 here means that this value will be named getLabel()
199 2 : readVesselKeywords(); // This makes sure resizing is done
200 2 : }
201 : }
202 :
203 : // Set the reference structure
204 62 : references.emplace_back( metricRegister().create<SingleDomainRMSD>( alignType ) );
205 62 : unsigned nn=references.size()-1;
206 62 : std::vector<double> align( structure.size(), 1.0 ), displace( structure.size(), 1.0 );
207 62 : references[nn]->setBoundsOnDistances( true, bondlength ); // We always use pbc
208 62 : references[nn]->setReferenceAtoms( structure, align, displace );
209 : // references[nn]->setNumberOfAtoms( structure.size() );
210 :
211 : // And prepare the task list
212 62 : deactivateAllTasks();
213 281340 : for(unsigned i=0; i<getFullNumberOfTasks(); ++i) {
214 281278 : taskFlags[i]=1;
215 : }
216 62 : lockContributors();
217 62 : }
218 :
219 2676 : void SecondaryStructureRMSD::calculate() {
220 2676 : runAllTasks();
221 2676 : }
222 :
223 929772 : void SecondaryStructureRMSD::performTask( const unsigned& task_index, const unsigned& current, MultiValue& myvals ) const {
224 : // Retrieve the positions
225 929772 : std::vector<Vector> pos( references[0]->getNumberOfAtoms() );
226 929772 : const unsigned n=pos.size();
227 28822932 : for(unsigned i=0; i<n; ++i) {
228 27893160 : pos[i]=ActionAtomistic::getPosition( getAtomIndex(current,i) );
229 : }
230 :
231 : // This does strands cutoff
232 929772 : Vector distance;
233 929772 : if( nopbc ) {
234 0 : distance=delta( pos[align_atom_1],pos[align_atom_2] );
235 : } else {
236 929772 : distance=pbcDistance( pos[align_atom_1],pos[align_atom_2] );
237 : }
238 929772 : if( s_cutoff2>0 ) {
239 927264 : if( distance.modulo2()>s_cutoff2 ) {
240 : myvals.setValue( 0, 0.0 );
241 842102 : return;
242 : }
243 : }
244 :
245 : // This aligns the two strands if this is required
246 87670 : if( alignType!="DRMSD" && align_strands && !nopbc ) {
247 928380 : for(unsigned i=0; i<14; ++i) {
248 866488 : const Vector & first (pos[i]);
249 866488 : Vector & second (pos[i+1]);
250 866488 : second=first+pbcDistance(first,second);
251 : }
252 866488 : for(unsigned i=16; i<n-1; ++i) {
253 804596 : const Vector & first (pos[i]);
254 804596 : Vector & second (pos[i+1]);
255 804596 : second=first+pbcDistance(first,second);
256 : }
257 61892 : Vector origin_old, origin_new;
258 61892 : origin_old=pos[align_atom_2];
259 61892 : origin_new=pos[align_atom_1]+distance;
260 990272 : for(unsigned i=15; i<30; ++i) {
261 928380 : pos[i]+=( origin_new - origin_old );
262 : }
263 25778 : } else if( alignType!="DRMSD" && !nopbc ) {
264 0 : for(unsigned i=0; i<n-1; ++i) {
265 0 : const Vector & first (pos[i]);
266 0 : Vector & second (pos[i+1]);
267 0 : second=first+pbcDistance(first,second);
268 : }
269 : }
270 : // Create a holder for the derivatives
271 87670 : ReferenceValuePack mypack( 0, pos.size(), myvals );
272 : mypack.setValIndex( 1 );
273 2717770 : for(unsigned i=0; i<n; ++i) {
274 : mypack.setAtomIndex( i, getAtomIndex(current,i) );
275 : }
276 :
277 : // And now calculate the RMSD
278 : const Pbc& pbc=getPbc();
279 : unsigned closest=0;
280 87670 : double r = references[0]->calculate( pos, pbc, mypack, false );
281 87670 : const unsigned rs = references.size();
282 130169 : for(unsigned i=1; i<rs; ++i) {
283 42499 : mypack.setValIndex( i+1 );
284 42499 : double nr=references[i]->calculate( pos, pbc, mypack, false );
285 42499 : if( nr<r ) {
286 : closest=i;
287 : r=nr;
288 : }
289 : }
290 :
291 : // Transfer everything to the value
292 : myvals.setValue( 0, 1.0 );
293 : myvals.setValue( 1, r );
294 87670 : if( closest>0 ) {
295 19914 : mypack.moveDerivatives( closest+1, 1 );
296 : }
297 :
298 87670 : if( !mypack.virialWasSet() ) {
299 61892 : Tensor vir;
300 61892 : const unsigned cacs = colvar_atoms[current].size();
301 1918652 : for(unsigned i=0; i<cacs; ++i) {
302 1856760 : vir+=(-1.0*Tensor( pos[i], mypack.getAtomDerivative(i) ));
303 : }
304 : mypack.setValIndex(1);
305 61892 : mypack.addBoxDerivatives( vir );
306 : }
307 :
308 : return;
309 87670 : }
310 :
311 300 : void SecondaryStructureRMSD::apply() {
312 300 : if( getForcesFromVessels( forcesToApply ) ) {
313 276 : setForcesOnAtoms( forcesToApply );
314 : }
315 300 : }
316 :
317 : }
318 : }
|
