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 28 : void SecondaryStructureRMSD::registerKeywords( Keywords& keys ) {
35 28 : Action::registerKeywords( keys );
36 28 : ActionWithValue::registerKeywords( keys );
37 28 : ActionAtomistic::registerKeywords( keys );
38 56 : 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 56 : 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 56 : keys.addFlag("NOPBC",false,"ignore the periodic boundary conditions");
49 56 : keys.add("compulsory","R_0","0.08","The r_0 parameter of the switching function.");
50 56 : keys.add("compulsory","D_0","0.0","The d_0 parameter of the switching function");
51 56 : keys.add("compulsory","NN","8","The n parameter of the switching function");
52 56 : keys.add("compulsory","MM","12","The m parameter of the switching function");
53 56 : 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 56 : keys.addFlag("VERBOSE",false,"write a more detailed output");
58 56 : 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 28 : ActionWithVessel::registerKeywords( keys );
61 140 : keys.use("LESS_THAN"); keys.use("MIN"); keys.use("ALT_MIN"); keys.use("LOWEST"); keys.use("HIGHEST");
62 28 : keys.setComponentsIntroduction("By default this Action calculates the number of structural units that are within a certain "
63 : "distance of a idealized secondary structure element. This quantity can then be referenced "
64 : "elsewhere in the input by using the label of the action. However, this Action can also be used to "
65 : "calculate the following quantities by using the keywords as described below. The quantities then "
66 : "calculated can be referenced using the label of the action followed by a dot and then the name "
67 : "from the table below. Please note that you can use the LESS_THAN keyword more than once. The resulting "
68 : "components will be labelled <em>label</em>.lessthan-1, <em>label</em>.lessthan-2 and so on unless you "
69 : "exploit the fact that these labels can be given custom labels by using the LABEL keyword in the "
70 : "description of you LESS_THAN function that you are computing");
71 28 : }
72 :
73 25 : SecondaryStructureRMSD::SecondaryStructureRMSD(const ActionOptions&ao):
74 : Action(ao),
75 : ActionAtomistic(ao),
76 : ActionWithValue(ao),
77 : ActionWithVessel(ao),
78 25 : nopbc(false),
79 25 : align_strands(false),
80 25 : s_cutoff2(0),
81 25 : align_atom_1(0),
82 25 : align_atom_2(0)
83 : {
84 50 : parse("TYPE",alignType); parseFlag("NOPBC",nopbc);
85 25 : log.printf(" distances from secondary structure elements are calculated using %s algorithm\n",alignType.c_str() );
86 75 : log<<" Bibliography "<<plumed.cite("Pietrucci and Laio, J. Chem. Theory Comput. 5, 2197 (2009)"); log<<"\n";
87 :
88 25 : parseFlag("VERBOSE",verbose_output);
89 :
90 50 : if( keywords.exists("STRANDS_CUTOFF") ) {
91 22 : double s_cutoff = 0;
92 22 : parse("STRANDS_CUTOFF",s_cutoff); align_strands=true;
93 22 : if( s_cutoff>0) log.printf(" ignoring contributions from strands that are more than %f apart\n",s_cutoff);
94 22 : s_cutoff2=s_cutoff*s_cutoff;
95 : }
96 25 : }
97 :
98 25 : SecondaryStructureRMSD::~SecondaryStructureRMSD() {
99 : // destructor needed to delete forward declarated objects
100 50 : }
101 :
102 62 : void SecondaryStructureRMSD::turnOnDerivatives() {
103 62 : ActionWithValue::turnOnDerivatives();
104 62 : needsDerivatives();
105 62 : }
106 :
107 22 : void SecondaryStructureRMSD::setAtomsFromStrands( const unsigned& atom1, const unsigned& atom2 ) {
108 22 : align_atom_1=atom1; align_atom_2=atom2;
109 22 : }
110 :
111 25 : void SecondaryStructureRMSD::readBackboneAtoms( const std::string& moltype, std::vector<unsigned>& chain_lengths ) {
112 25 : auto* moldat=plumed.getActionSet().selectLatest<GenericMolInfo*>(this);
113 25 : if( ! moldat ) error("Unable to find MOLINFO in input");
114 :
115 25 : std::vector<std::string> resstrings; parseVector( "RESIDUES", resstrings );
116 25 : if( !verbose_output ) {
117 25 : if(resstrings.size()==0) error("residues are not defined, check the keyword RESIDUES");
118 25 : else if(resstrings[0]=="all") {
119 21 : log.printf(" examining all possible secondary structure combinations\n");
120 : } else {
121 4 : log.printf(" examining secondary structure in residue positions : %s \n",resstrings[0].c_str() );
122 6 : for(unsigned i=1; i<resstrings.size(); ++i) log.printf(", %s",resstrings[i].c_str() );
123 4 : log.printf("\n");
124 : }
125 : }
126 : std::vector< std::vector<AtomNumber> > backatoms;
127 25 : moldat->getBackbone( resstrings, moltype, backatoms );
128 :
129 25 : chain_lengths.resize( backatoms.size() );
130 358 : for(unsigned i=0; i<backatoms.size(); ++i) {
131 333 : chain_lengths[i]=backatoms[i].size();
132 13593 : for(unsigned j=0; j<backatoms[i].size(); ++j) all_atoms.push_back( backatoms[i][j] );
133 : }
134 25 : ActionAtomistic::requestAtoms( all_atoms );
135 25 : forcesToApply.resize( getNumberOfDerivatives() );
136 25 : }
137 :
138 99342 : void SecondaryStructureRMSD::addColvar( const std::vector<unsigned>& newatoms ) {
139 99342 : if( colvar_atoms.size()>0 ) plumed_assert( colvar_atoms[0].size()==newatoms.size() );
140 99342 : if( verbose_output ) {
141 0 : log.printf(" Secondary structure segment %u contains atoms : ", static_cast<unsigned>(colvar_atoms.size()+1));
142 0 : for(unsigned i=0; i<newatoms.size(); ++i) log.printf("%d ",all_atoms[newatoms[i]].serial() );
143 0 : log.printf("\n");
144 : }
145 99342 : addTaskToList( colvar_atoms.size() );
146 99342 : colvar_atoms.push_back( newatoms );
147 99342 : }
148 :
149 35 : void SecondaryStructureRMSD::setSecondaryStructure( std::vector<Vector>& structure, double bondlength, double units ) {
150 : // If we are in natural units get conversion factor from nm into natural length units
151 35 : if( plumed.getAtoms().usingNaturalUnits() ) {
152 0 : error("cannot use this collective variable when using natural units");
153 : }
154 35 : plumed_massert( !(align_strands && align_atom_1==0 && align_atom_2==0), "you must use setAtomsFromStrands with strands cutoff");
155 :
156 : // Convert into correct units
157 1085 : for(unsigned i=0; i<structure.size(); ++i) {
158 1050 : structure[i][0]*=units; structure[i][1]*=units; structure[i][2]*=units;
159 : }
160 :
161 35 : if( references.size()==0 ) {
162 25 : readVesselKeywords();
163 25 : if( getNumberOfVessels()==0 ) {
164 4 : double r0; parse("R_0",r0); double d0; parse("D_0",d0);
165 4 : int nn; parse("NN",nn); int mm; parse("MM",mm);
166 2 : std::ostringstream ostr;
167 2 : ostr<<"RATIONAL R_0="<<r0<<" D_0="<<d0<<" NN="<<nn<<" MM="<<mm;
168 2 : std::string input=ostr.str(); addVessel( "LESS_THAN", input, -1 ); // -1 here means that this value will be named getLabel()
169 2 : readVesselKeywords(); // This makes sure resizing is done
170 2 : }
171 : }
172 :
173 : // Set the reference structure
174 35 : references.emplace_back( metricRegister().create<SingleDomainRMSD>( alignType ) );
175 35 : unsigned nn=references.size()-1;
176 35 : std::vector<double> align( structure.size(), 1.0 ), displace( structure.size(), 1.0 );
177 35 : references[nn]->setBoundsOnDistances( true, bondlength ); // We always use pbc
178 35 : references[nn]->setReferenceAtoms( structure, align, displace );
179 : // references[nn]->setNumberOfAtoms( structure.size() );
180 :
181 : // And prepare the task list
182 35 : deactivateAllTasks();
183 148959 : for(unsigned i=0; i<getFullNumberOfTasks(); ++i) taskFlags[i]=1;
184 35 : lockContributors();
185 35 : }
186 :
187 2568 : void SecondaryStructureRMSD::calculate() {
188 2568 : runAllTasks();
189 2568 : }
190 :
191 400032 : void SecondaryStructureRMSD::performTask( const unsigned& task_index, const unsigned& current, MultiValue& myvals ) const {
192 : // Retrieve the positions
193 400032 : std::vector<Vector> pos( references[0]->getNumberOfAtoms() );
194 400032 : const unsigned n=pos.size();
195 12400992 : for(unsigned i=0; i<n; ++i) pos[i]=ActionAtomistic::getPosition( getAtomIndex(current,i) );
196 :
197 : // This does strands cutoff
198 400032 : Vector distance;
199 400032 : if( nopbc ) distance=delta( pos[align_atom_1],pos[align_atom_2] );
200 400032 : else distance=pbcDistance( pos[align_atom_1],pos[align_atom_2] );
201 400032 : if( s_cutoff2>0 ) {
202 397524 : if( distance.modulo2()>s_cutoff2 ) {
203 : myvals.setValue( 0, 0.0 );
204 360914 : return;
205 : }
206 : }
207 :
208 : // This aligns the two strands if this is required
209 39118 : if( alignType!="DRMSD" && align_strands && !nopbc ) {
210 382260 : for(unsigned i=0; i<14; ++i) {
211 356776 : const Vector & first (pos[i]);
212 356776 : Vector & second (pos[i+1]);
213 356776 : second=first+pbcDistance(first,second);
214 : }
215 356776 : for(unsigned i=16; i<n-1; ++i) {
216 331292 : const Vector & first (pos[i]);
217 331292 : Vector & second (pos[i+1]);
218 331292 : second=first+pbcDistance(first,second);
219 : }
220 25484 : Vector origin_old, origin_new; origin_old=pos[align_atom_2];
221 25484 : origin_new=pos[align_atom_1]+distance;
222 407744 : for(unsigned i=15; i<30; ++i) {
223 382260 : pos[i]+=( origin_new - origin_old );
224 : }
225 13634 : } else if( alignType!="DRMSD" && !nopbc ) {
226 0 : for(unsigned i=0; i<n-1; ++i) {
227 0 : const Vector & first (pos[i]);
228 0 : Vector & second (pos[i+1]);
229 0 : second=first+pbcDistance(first,second);
230 : }
231 : }
232 : // Create a holder for the derivatives
233 39118 : ReferenceValuePack mypack( 0, pos.size(), myvals ); mypack.setValIndex( 1 );
234 1212658 : for(unsigned i=0; i<n; ++i) mypack.setAtomIndex( i, getAtomIndex(current,i) );
235 :
236 : // And now calculate the RMSD
237 : const Pbc& pbc=getPbc();
238 : unsigned closest=0;
239 39118 : double r = references[0]->calculate( pos, pbc, mypack, false );
240 39118 : const unsigned rs = references.size();
241 57377 : for(unsigned i=1; i<rs; ++i) {
242 18259 : mypack.setValIndex( i+1 );
243 18259 : double nr=references[i]->calculate( pos, pbc, mypack, false );
244 18259 : if( nr<r ) { closest=i; r=nr; }
245 : }
246 :
247 : // Transfer everything to the value
248 : myvals.setValue( 0, 1.0 ); myvals.setValue( 1, r );
249 39118 : if( closest>0 ) mypack.moveDerivatives( closest+1, 1 );
250 :
251 39118 : if( !mypack.virialWasSet() ) {
252 25484 : Tensor vir;
253 25484 : const unsigned cacs = colvar_atoms[current].size();
254 790004 : for(unsigned i=0; i<cacs; ++i) {
255 764520 : vir+=(-1.0*Tensor( pos[i], mypack.getAtomDerivative(i) ));
256 : }
257 25484 : mypack.setValIndex(1); mypack.addBoxDerivatives( vir );
258 : }
259 :
260 : return;
261 39118 : }
262 :
263 192 : void SecondaryStructureRMSD::apply() {
264 192 : if( getForcesFromVessels( forcesToApply ) ) setForcesOnAtoms( forcesToApply );
265 192 : }
266 :
267 : }
268 : }
|
