Line data Source code
1 : /* +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
2 : Copyright (c) 2013-2019 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 112 : 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 140 : 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 140 : keys.add("compulsory","R_0","0.08","The r_0 parameter of the switching function.");
49 140 : keys.add("compulsory","D_0","0.0","The d_0 parameter of the switching function");
50 140 : keys.add("compulsory","NN","8","The n parameter of the switching function");
51 140 : keys.add("compulsory","MM","12","The m parameter of the switching function");
52 112 : 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 : "However, if you are using some other option, then this cannot be used");
56 84 : keys.addFlag("VERBOSE",false,"write a more detailed output");
57 112 : keys.add("hidden","NL_STRIDE","the frequency with which the neighbor list should be updated. Between neighbour list update steps all quantities "
58 : "that contributed less than TOL at the previous neighbor list update step are ignored.");
59 28 : ActionWithVessel::registerKeywords( keys );
60 168 : keys.use("LESS_THAN"); keys.use("MIN"); keys.use("ALT_MIN"); keys.use("LOWEST"); keys.use("HIGHEST");
61 56 : 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 : "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 75 : align_atom_2(0)
81 : {
82 50 : parse("TYPE",alignType);
83 50 : log.printf(" distances from secondary structure elements are calculated using %s algorithm\n",alignType.c_str() );
84 75 : log<<" Bibliography "<<plumed.cite("Pietrucci and Laio, J. Chem. Theory Comput. 5, 2197 (2009)"); log<<"\n";
85 :
86 50 : parseFlag("VERBOSE",verbose_output);
87 :
88 50 : if( keywords.exists("STRANDS_CUTOFF") ) {
89 22 : double s_cutoff = 0;
90 44 : 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 75 : SecondaryStructureRMSD::~SecondaryStructureRMSD() {
97 155 : for(unsigned i=0; i<references.size(); ++i) delete references[i];
98 25 : }
99 :
100 62 : void SecondaryStructureRMSD::turnOnDerivatives() {
101 62 : ActionWithValue::turnOnDerivatives();
102 62 : needsDerivatives();
103 62 : }
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 50 : std::vector<SetupMolInfo*> moldat=plumed.getActionSet().select<SetupMolInfo*>();
111 25 : if( moldat.size()==0 ) error("Unable to find MOLINFO in input");
112 :
113 75 : std::vector<std::string> resstrings; parseVector( "RESIDUES", resstrings );
114 25 : if( !verbose_output ) {
115 25 : if(resstrings.size()==0) error("residues are not defined, check the keyword RESIDUES");
116 25 : else if(resstrings[0]=="all") {
117 21 : log.printf(" examining all possible secondary structure combinations\n");
118 : } else {
119 8 : log.printf(" examining secondary structure in residue positions : %s \n",resstrings[0].c_str() );
120 14 : for(unsigned i=1; i<resstrings.size(); ++i) log.printf(", %s",resstrings[i].c_str() );
121 4 : log.printf("\n");
122 : }
123 : }
124 25 : std::vector< std::vector<AtomNumber> > backatoms;
125 25 : moldat[0]->getBackbone( resstrings, moltype, backatoms );
126 :
127 25 : chain_lengths.resize( backatoms.size() );
128 1049 : for(unsigned i=0; i<backatoms.size(); ++i) {
129 333 : chain_lengths[i]=backatoms[i].size();
130 40446 : for(unsigned j=0; j<backatoms[i].size(); ++j) all_atoms.push_back( backatoms[i][j] );
131 : }
132 25 : ActionAtomistic::requestAtoms( all_atoms );
133 25 : forcesToApply.resize( getNumberOfDerivatives() );
134 25 : }
135 :
136 99342 : void SecondaryStructureRMSD::addColvar( const std::vector<unsigned>& newatoms ) {
137 198659 : if( colvar_atoms.size()>0 ) plumed_assert( colvar_atoms[0].size()==newatoms.size() );
138 99342 : if( verbose_output ) {
139 0 : log.printf(" Secondary structure segment %u contains atoms : ", static_cast<unsigned>(colvar_atoms.size()+1));
140 0 : for(unsigned i=0; i<newatoms.size(); ++i) log.printf("%d ",all_atoms[newatoms[i]].serial() );
141 0 : log.printf("\n");
142 : }
143 198684 : addTaskToList( colvar_atoms.size() );
144 99342 : colvar_atoms.push_back( newatoms );
145 99342 : }
146 :
147 35 : void SecondaryStructureRMSD::setSecondaryStructure( std::vector<Vector>& structure, double bondlength, double units ) {
148 : // If we are in natural units get conversion factor from nm into natural length units
149 35 : if( plumed.getAtoms().usingNaturalUnits() ) {
150 0 : error("cannot use this collective variable when using natural units");
151 : }
152 35 : plumed_massert( !(align_strands && align_atom_1==0 && align_atom_2==0), "you must use setAtomsFromStrands with strands cutoff");
153 :
154 : // Convert into correct units
155 3220 : for(unsigned i=0; i<structure.size(); ++i) {
156 3150 : structure[i][0]*=units; structure[i][1]*=units; structure[i][2]*=units;
157 : }
158 :
159 35 : if( references.size()==0 ) {
160 25 : readVesselKeywords();
161 25 : if( getNumberOfVessels()==0 ) {
162 6 : double r0; parse("R_0",r0); double d0; parse("D_0",d0);
163 6 : int nn; parse("NN",nn); int mm; parse("MM",mm);
164 4 : std::ostringstream ostr;
165 8 : ostr<<"RATIONAL R_0="<<r0<<" D_0="<<d0<<" NN="<<nn<<" MM="<<mm;
166 4 : std::string input=ostr.str(); addVessel( "LESS_THAN", input, -1 ); // -1 here means that this value will be named getLabel()
167 2 : readVesselKeywords(); // This makes sure resizing is done
168 : }
169 : }
170 :
171 : // Set the reference structure
172 70 : references.push_back( metricRegister().create<SingleDomainRMSD>( alignType ) );
173 35 : unsigned nn=references.size()-1;
174 105 : std::vector<double> align( structure.size(), 1.0 ), displace( structure.size(), 1.0 );
175 70 : references[nn]->setBoundsOnDistances( true, bondlength ); // We always use pbc
176 35 : references[nn]->setReferenceAtoms( structure, align, displace );
177 : // references[nn]->setNumberOfAtoms( structure.size() );
178 :
179 : // And prepare the task list
180 35 : deactivateAllTasks();
181 297883 : for(unsigned i=0; i<getFullNumberOfTasks(); ++i) taskFlags[i]=1;
182 35 : lockContributors();
183 35 : }
184 :
185 2568 : void SecondaryStructureRMSD::calculate() {
186 2568 : runAllTasks();
187 2568 : }
188 :
189 400032 : void SecondaryStructureRMSD::performTask( const unsigned& task_index, const unsigned& current, MultiValue& myvals ) const {
190 : // Retrieve the positions
191 800064 : std::vector<Vector> pos( references[0]->getNumberOfAtoms() );
192 400032 : const unsigned n=pos.size();
193 36402912 : for(unsigned i=0; i<n; ++i) pos[i]=ActionAtomistic::getPosition( getAtomIndex(current,i) );
194 :
195 : // This does strands cutoff
196 1200096 : Vector distance=pbcDistance( pos[align_atom_1],pos[align_atom_2] );
197 400032 : if( s_cutoff2>0 ) {
198 397524 : if( distance.modulo2()>s_cutoff2 ) {
199 : myvals.setValue( 0, 0.0 );
200 360914 : return;
201 : }
202 : }
203 :
204 : // This aligns the two strands if this is required
205 78236 : if( alignType!="DRMSD" && align_strands ) {
206 50968 : Vector origin_old, origin_new; origin_old=pos[align_atom_2];
207 50968 : origin_new=pos[align_atom_1]+distance;
208 790004 : for(unsigned i=15; i<30; ++i) {
209 764520 : pos[i]+=( origin_new - origin_old );
210 : }
211 : }
212 : // Create a holder for the derivatives
213 78236 : ReferenceValuePack mypack( 0, pos.size(), myvals ); mypack.setValIndex( 1 );
214 2386198 : for(unsigned i=0; i<n; ++i) mypack.setAtomIndex( i, getAtomIndex(current,i) );
215 :
216 : // And now calculate the RMSD
217 : const Pbc& pbc=getPbc();
218 : unsigned closest=0;
219 39118 : double r = references[0]->calculate( pos, pbc, mypack, false );
220 39118 : const unsigned rs = references.size();
221 57377 : for(unsigned i=1; i<rs; ++i) {
222 18259 : mypack.setValIndex( i+1 );
223 36518 : double nr=references[i]->calculate( pos, pbc, mypack, false );
224 18259 : if( nr<r ) { closest=i; r=nr; }
225 : }
226 :
227 : // Transfer everything to the value
228 : myvals.setValue( 0, 1.0 ); myvals.setValue( 1, r );
229 39118 : if( closest>0 ) mypack.moveDerivatives( closest+1, 1 );
230 :
231 39118 : if( !mypack.virialWasSet() ) {
232 25484 : Tensor vir;
233 50968 : const unsigned cacs = colvar_atoms[current].size();
234 790004 : for(unsigned i=0; i<cacs; ++i) {
235 1529040 : vir+=(-1.0*Tensor( pos[i], mypack.getAtomDerivative(i) ));
236 : }
237 25484 : mypack.setValIndex(1); mypack.addBoxDerivatives( vir );
238 : }
239 :
240 : return;
241 : }
242 :
243 192 : void SecondaryStructureRMSD::apply() {
244 192 : if( getForcesFromVessels( forcesToApply ) ) setForcesOnAtoms( forcesToApply );
245 192 : }
246 :
247 : }
248 4839 : }
|
