Line data Source code
1 : /* +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
2 : Copyright (c) 2013-2018 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/SetupMolInfo.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 28 : void SecondaryStructureRMSD::registerKeywords( Keywords& keys ) {
35 28 : Action::registerKeywords( keys );
36 28 : ActionWithValue::registerKeywords( keys );
37 28 : ActionAtomistic::registerKeywords( keys );
38 : 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 28 : "As such if you define portions of the chain with fewer than N residues the code will crash.");
45 : 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 28 : "DRMSD method see \\ref DRMSD.");
48 28 : keys.add("compulsory","R_0","The r_0 parameter of the switching function.");
49 28 : keys.add("compulsory","D_0","0.0","The d_0 parameter of the switching function");
50 28 : keys.add("compulsory","NN","8","The n parameter of the switching function");
51 28 : keys.add("compulsory","MM","12","The m parameter of the switching function");
52 : keys.reserve("optional","STRANDS_CUTOFF","If in a segment of protein the two strands are further apart then the calculation "
53 : "of the actual RMSD is skipped as the structure is very far from being beta-sheet like. "
54 : "This keyword speeds up the calculation enormously when you are using the LESS_THAN option. "
55 28 : "However, if you are using some other option, then this cannot be used");
56 28 : keys.addFlag("VERBOSE",false,"write a more detailed output");
57 : keys.add("hidden","NL_STRIDE","the frequency with which the neighbor list should be updated. Between neighbour list update steps all quantities "
58 28 : "that contributed less than TOL at the previous neighbor list update step are ignored.");
59 28 : ActionWithVessel::registerKeywords( keys );
60 28 : keys.use("LESS_THAN"); keys.use("MIN"); keys.use("ALT_MIN"); keys.use("LOWEST"); keys.use("HIGHEST");
61 : keys.setComponentsIntroduction("By default this Action calculates the number of structural units that are within a certain "
62 : "distance of a idealised secondary structure element. This quantity can then be referenced "
63 : "elsewhere in the input by using the label of the action. However, this Action can also be used to "
64 : "calculate the following quantities by using the keywords as described below. The quantities then "
65 : "calculated can be referened using the label of the action followed by a dot and then the name "
66 : "from the table below. Please note that you can use the LESS_THAN keyword more than once. The resulting "
67 : "components will be labelled <em>label</em>.lessthan-1, <em>label</em>.lessthan-2 and so on unless you "
68 : "exploit the fact that these labels are customizable. In particular, by using the LABEL keyword in the "
69 28 : "description of you LESS_THAN function you can set name of the component that you are calculating");
70 28 : }
71 :
72 25 : SecondaryStructureRMSD::SecondaryStructureRMSD(const ActionOptions&ao):
73 : Action(ao),
74 : ActionAtomistic(ao),
75 : ActionWithValue(ao),
76 : ActionWithVessel(ao),
77 : align_strands(false),
78 : s_cutoff2(0),
79 : align_atom_1(0),
80 25 : align_atom_2(0)
81 : {
82 25 : parse("TYPE",alignType);
83 25 : log.printf(" distances from secondary structure elements are calculated using %s algorithm\n",alignType.c_str() );
84 25 : log<<" Bibliography "<<plumed.cite("Pietrucci and Laio, J. Chem. Theory Comput. 5, 2197 (2009)"); log<<"\n";
85 :
86 25 : parseFlag("VERBOSE",verbose_output);
87 :
88 25 : if( keywords.exists("STRANDS_CUTOFF") ) {
89 22 : double s_cutoff = 0;
90 22 : parse("STRANDS_CUTOFF",s_cutoff); align_strands=true;
91 22 : if( s_cutoff>0) log.printf(" ignoring contributions from strands that are more than %f apart\n",s_cutoff);
92 22 : s_cutoff2=s_cutoff*s_cutoff;
93 : }
94 25 : }
95 :
96 50 : SecondaryStructureRMSD::~SecondaryStructureRMSD() {
97 25 : for(unsigned i=0; i<references.size(); ++i) delete references[i];
98 25 : }
99 :
100 52 : void SecondaryStructureRMSD::turnOnDerivatives() {
101 52 : ActionWithValue::turnOnDerivatives();
102 52 : needsDerivatives();
103 52 : }
104 :
105 22 : void SecondaryStructureRMSD::setAtomsFromStrands( const unsigned& atom1, const unsigned& atom2 ) {
106 22 : align_atom_1=atom1; align_atom_2=atom2;
107 22 : }
108 :
109 25 : void SecondaryStructureRMSD::readBackboneAtoms( const std::string& moltype, std::vector<unsigned>& chain_lengths ) {
110 25 : std::vector<SetupMolInfo*> moldat=plumed.getActionSet().select<SetupMolInfo*>();
111 25 : if( moldat.size()==0 ) error("Unable to find MOLINFO in input");
112 :
113 50 : std::vector<std::string> resstrings; parseVector( "RESIDUES", resstrings );
114 25 : if( !verbose_output ) {
115 25 : if(resstrings[0]=="all") {
116 21 : log.printf(" examining all possible secondary structure combinations\n");
117 : } else {
118 4 : log.printf(" examining secondary structure in residue positions : %s \n",resstrings[0].c_str() );
119 4 : for(unsigned i=1; i<resstrings.size(); ++i) log.printf(", %s",resstrings[i].c_str() );
120 4 : log.printf("\n");
121 : }
122 : }
123 50 : std::vector< std::vector<AtomNumber> > backatoms;
124 25 : moldat[0]->getBackbone( resstrings, moltype, backatoms );
125 :
126 25 : chain_lengths.resize( backatoms.size() );
127 358 : for(unsigned i=0; i<backatoms.size(); ++i) {
128 333 : chain_lengths[i]=backatoms[i].size();
129 333 : for(unsigned j=0; j<backatoms[i].size(); ++j) all_atoms.push_back( backatoms[i][j] );
130 : }
131 25 : ActionAtomistic::requestAtoms( all_atoms );
132 50 : forcesToApply.resize( getNumberOfDerivatives() );
133 25 : }
134 :
135 99342 : void SecondaryStructureRMSD::addColvar( const std::vector<unsigned>& newatoms ) {
136 99342 : if( colvar_atoms.size()>0 ) plumed_assert( colvar_atoms[0].size()==newatoms.size() );
137 99342 : if( verbose_output ) {
138 0 : log.printf(" Secondary structure segment %u contains atoms : ", static_cast<unsigned>(colvar_atoms.size()+1));
139 0 : for(unsigned i=0; i<newatoms.size(); ++i) log.printf("%d ",all_atoms[newatoms[i]].serial() );
140 0 : log.printf("\n");
141 : }
142 99342 : addTaskToList( colvar_atoms.size() );
143 99342 : colvar_atoms.push_back( newatoms );
144 99342 : }
145 :
146 35 : void SecondaryStructureRMSD::setSecondaryStructure( std::vector<Vector>& structure, double bondlength, double units ) {
147 : // If we are in natural units get conversion factor from nm into natural length units
148 35 : if( plumed.getAtoms().usingNaturalUnits() ) {
149 0 : error("cannot use this collective variable when using natural units");
150 : }
151 35 : plumed_massert( !(align_strands && align_atom_1==0 && align_atom_2==0), "you must use setAtomsFromStrands with strands cutoff");
152 :
153 : // Convert into correct units
154 1085 : for(unsigned i=0; i<structure.size(); ++i) {
155 1050 : structure[i][0]*=units; structure[i][1]*=units; structure[i][2]*=units;
156 : }
157 :
158 35 : if( references.size()==0 ) {
159 25 : readVesselKeywords();
160 25 : if( getNumberOfVessels()==0 ) {
161 2 : double r0; parse("R_0",r0); double d0; parse("D_0",d0);
162 2 : int nn; parse("NN",nn); int mm; parse("MM",mm);
163 2 : std::ostringstream ostr;
164 2 : ostr<<"RATIONAL R_0="<<r0<<" D_0="<<d0<<" NN="<<nn<<" MM="<<mm;
165 4 : std::string input=ostr.str(); addVessel( "LESS_THAN", input, -1 ); // -1 here means that this value will be named getLabel()
166 4 : readVesselKeywords(); // This makes sure resizing is done
167 : }
168 : }
169 :
170 : // Set the reference structure
171 35 : references.push_back( metricRegister().create<SingleDomainRMSD>( alignType ) );
172 35 : unsigned nn=references.size()-1;
173 70 : std::vector<double> align( structure.size(), 1.0 ), displace( structure.size(), 1.0 );
174 35 : references[nn]->setBoundsOnDistances( true, bondlength ); // We always use pbc
175 35 : references[nn]->setReferenceAtoms( structure, align, displace );
176 : // references[nn]->setNumberOfAtoms( structure.size() );
177 :
178 : // And prepare the task list
179 35 : deactivateAllTasks();
180 35 : for(unsigned i=0; i<getFullNumberOfTasks(); ++i) taskFlags[i]=1;
181 70 : lockContributors();
182 35 : }
183 :
184 2568 : void SecondaryStructureRMSD::calculate() {
185 2568 : runAllTasks();
186 2568 : }
187 :
188 400018 : void SecondaryStructureRMSD::performTask( const unsigned& task_index, const unsigned& current, MultiValue& myvals ) const {
189 : // Retrieve the positions
190 400018 : std::vector<Vector> pos( references[0]->getNumberOfAtoms() );
191 400029 : const unsigned n=pos.size();
192 397741 : for(unsigned i=0; i<n; ++i) pos[i]=ActionAtomistic::getPosition( getAtomIndex(current,i) );
193 :
194 : // This does strands cutoff
195 399983 : Vector distance=pbcDistance( pos[align_atom_1],pos[align_atom_2] );
196 400027 : if( s_cutoff2>0 ) {
197 397518 : if( distance.modulo2()>s_cutoff2 ) {
198 360905 : myvals.setValue( 0, 0.0 );
199 360912 : return;
200 : }
201 : }
202 :
203 : // This aligns the two strands if this is required
204 39092 : if( alignType!="DRMSD" && align_strands ) {
205 25482 : Vector origin_old, origin_new; origin_old=pos[align_atom_2];
206 25482 : origin_new=pos[align_atom_1]+distance;
207 407611 : for(unsigned i=15; i<30; ++i) {
208 382128 : pos[i]+=( origin_new - origin_old );
209 : }
210 : }
211 : // Create a holder for the derivatives
212 78229 : ReferenceValuePack mypack( 0, pos.size(), myvals ); mypack.setValIndex( 1 );
213 38884 : for(unsigned i=0; i<n; ++i) mypack.setAtomIndex( i, getAtomIndex(current,i) );
214 :
215 : // And now calculate the RMSD
216 39115 : const Pbc& pbc=getPbc();
217 39109 : unsigned closest=0;
218 39109 : double r = references[0]->calculate( pos, pbc, mypack, false );
219 39109 : const unsigned rs = references.size();
220 57368 : for(unsigned i=1; i<rs; ++i) {
221 18253 : mypack.setValIndex( i+1 );
222 18252 : double nr=references[i]->calculate( pos, pbc, mypack, false );
223 18258 : if( nr<r ) { closest=i; r=nr; }
224 : }
225 :
226 : // Transfer everything to the value
227 39115 : myvals.setValue( 0, 1.0 ); myvals.setValue( 1, r );
228 39114 : if( closest>0 ) mypack.moveDerivatives( closest+1, 1 );
229 :
230 39113 : if( !mypack.virialWasSet() ) {
231 25481 : Tensor vir;
232 25484 : const unsigned cacs = colvar_atoms[current].size();
233 789932 : for(unsigned i=0; i<cacs; ++i) {
234 764498 : vir+=(-1.0*Tensor( pos[i], mypack.getAtomDerivative(i) ));
235 : }
236 25434 : mypack.setValIndex(1); mypack.addBoxDerivatives( vir );
237 : }
238 :
239 439138 : return;
240 : }
241 :
242 192 : void SecondaryStructureRMSD::apply() {
243 192 : if( getForcesFromVessels( forcesToApply ) ) setForcesOnAtoms( forcesToApply );
244 192 : }
245 :
246 : }
247 2523 : }
|
