Line data Source code
1 : /* +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
2 : Copyright (c) 2017-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 "SecondaryStructureBase.h"
23 : #include "core/ActionRegister.h"
24 : #include "tools/RMSD.h"
25 :
26 : //+PLUMEDOC MCOLVAR SECONDARY_STRUCTURE_RMSD
27 : /*
28 : Calclulate the RMSD between segments of a protein and a reference structure of interest
29 :
30 : This action is used in the shortcuts [ALPHARMSD](ALPHARMSD.md), [ANTIBETARMSD](ANTIBETARMSD.md) and [PARABETARMSD](PARABETARMSD.md). It calculates a
31 : vector of [RMSD](RMSD.md) values between a single reference multiple configurations and the instantaneous
32 : positions of various groups of atoms. For example, in the following input we define a single set of reference
33 : set of coordinates for 3 atoms.
34 :
35 : ```plumed
36 : c1: SECONDARY_STRUCTURE_RMSD ...
37 : ATOMS=13-24
38 : STRUCTURE1=1,0,0,0,1,0,0,0,1
39 : SEGMENT1=0,1,2 SEGMENT2=3,4,5
40 : SEGMENT3=6,7,8 SEGMENT4=9,10,11
41 : TYPE=OPTIMAL
42 : ...
43 : PRINT ARG=c1 FILE=colvar
44 : ```
45 :
46 : A four dimensional vector is then returned that contains the RMSD distances between the 4 sets of 3 atoms that were specified using the `SEGMENT` keywords
47 : and the reference coordinates. Notice that you can use multiple instances of the `STRUCTURE` keyword. In general the number of vectors output
48 : is equal to the number of times the `STRUCTURE` keyword is used. Further note that the indices specified to the SEGMENT keywords give the indices in the list
49 : of atoms specified by the ATOMS keyword. `SEGMENT1=0,1,2` in the above input is measuring the distance between the instantaneous positions of atoms 13, 14 and 15
50 : and the reference structure.
51 :
52 : ## Periodic boundary conditions
53 :
54 : If you turn off periodic boundary conditions by using the NOPBC flag as shown below:
55 :
56 : ```plumed
57 : c1: SECONDARY_STRUCTURE_RMSD ...
58 : ATOMS=13-24 NOPBC
59 : STRUCTURE1=1,0,0,0,1,0,0,0,1
60 : SEGMENT1=0,1,2 SEGMENT2=3,4,5
61 : SEGMENT3=6,7,8 SEGMENT4=9,10,11
62 : TYPE=SIMPLE
63 : ...
64 : PRINT ARG=c1 FILE=colvar
65 : ```
66 :
67 : the distances between in the instantaneous structure are evaluated in a way that does not take the periodic boundary conditions
68 : into account. Whenver this flag is __not__ present the distances in the instantaneous structure are evaluated in a way that takes the periodic boundary conditions
69 : into account.
70 :
71 : If you use the `ALIGN_STRANDS` flag to evaluate the distances in a way that takes the periodic boundary conditions into accounts means that
72 : calculating the distance between the positions of the 7th and 22nd atoms ($\mathbf{x}_7$ and $\mathbf{x}_{22}$) of the instanenous structures
73 : using the periodic boundary conditions. The 22 atom is the repositioned at:
74 :
75 : $$
76 : \mathbf{x}_{22} = \mathbf{x}_7 + ( \mathbf{x}_{22} - \mathbf{x}_7 )
77 : $$
78 :
79 : where the difference in the second term on the right of the equality sign is evaluated using the periodic boundary conditions. New positions are then determined
80 : for atoms 16 through 30 by adding the difference (evaluated without PBC) between the new position for atom 22 after the above transformation has been performed and the
81 : original position of this atom. This treatment is designed for making sure that the two strands of the [ANTIBETARMSD](ANTIBETARMSD.md) and [PARABETARMSD](PARABETARMSD.md)
82 : sheets are not broken by the periodic boundaries.
83 :
84 : ## Working with beta sheet like structures
85 :
86 : As discussed in the documentation for [ANTIBETARMSD](ANTIBETARMSD.md) and [PARABETARMSD](PARABETARMSD.md) using the STRANDS_KEYWORD cutoff can speed up the calculation
87 : dramatically. The reason this keyword gives such dramatric speed ups is illustrated by the example input:
88 :
89 : ```plumed
90 : ab_cut_dists: DISTANCE ...
91 : ATOMS1=19,69 ATOMS2=19,79
92 : ATOMS3=19,89 ATOMS4=19,99
93 : ATOMS5=19,109
94 : ...
95 : ab_cut: CUSTOM ARG=ab_cut_dists FUNC=step(1-x) PERIODIC=NO
96 : ab_rmsd: SECONDARY_STRUCTURE_RMSD ...
97 : TYPE=OPTIMAL ALIGN_STRANDS MASK=ab_cut
98 : ATOMS=7,9,11,15,16,17,19,21,25,26,27,29,31,35,36,37,39,41,45,46,47,49,51,55,56,57,59,61,65,66,67,69,71,75,76,77,79,81,85,86,87,89,91,95,96,97,99,101,105,106,107,109,111,115,116,117,119,121,125,126
99 : SEGMENT1=0,1,2,3,4,5,6,7,8,9,10,11,12,13,14,25,26,27,28,29,30,31,32,33,34,35,36,37,38,39
100 : SEGMENT2=0,1,2,3,4,5,6,7,8,9,10,11,12,13,14,30,31,32,33,34,35,36,37,38,39,40,41,42,43,44
101 : SEGMENT3=0,1,2,3,4,5,6,7,8,9,10,11,12,13,14,35,36,37,38,39,40,41,42,43,44,45,46,47,48,49
102 : SEGMENT4=0,1,2,3,4,5,6,7,8,9,10,11,12,13,14,40,41,42,43,44,45,46,47,48,49,50,51,52,53,54
103 : SEGMENT5=0,1,2,3,4,5,6,7,8,9,10,11,12,13,14,45,46,47,48,49,50,51,52,53,54,55,56,57,58,59
104 : STRUCTURE1=2.263,-3.795,1.722,2.493,-2.426,2.263,3.847,-1.838,1.761,1.301,-1.517,1.921,0.852,-1.504,0.739,0.818,-0.738,2.917,-0.299,0.243,2.748,-1.421,-0.076,3.757,0.273,1.68,2.854,0.902,1.993,3.888,0.119,2.532,1.813,0.683,3.916,1.68,1.58,3.94,0.395,-0.394,5.011,1.63,-1.459,4.814,0.982,-2.962,3.559,-1.359,-2.439,2.526,-2.287,-1.189,3.006,-3.087,-2.081,1.231,-1.52,-1.524,1.324,-0.409,-2.326,0.037,-2.095,-1.858,-1.269,-1.554,-3.053,-2.199,-1.291,-0.869,-1.949,-2.512,-1.255,-2.07,-3.71,0.326,-2.363,-2.072,1.405,-2.992,-2.872,2.699,-2.129,-2.917,1.745,-4.399,-2.33,1.899,-4.545,-1.102
105 : ...
106 : ab_ltu: LESS_THAN ARG=ab_rmsd SWITCH={RATIONAL R_0=0.1 D_0=0.0 NN=8 MM=12} MASK=ab_cut
107 : ab_lt: CUSTOM ARG=ab_ltu,ab_cut FUNC=x*y PERIODIC=NO
108 : ab: SUM ARG=ab_lt PERIODIC=NO
109 : PRINT ARG=ab FILE=colvar
110 : ```
111 :
112 : The input above illustrates the input that is used by the shortcut that computes the [ANTIBETARMSD](ANTIBETARMSD.md). You can see that the distances between the atoms on the two strands of the
113 : beta sheet are computed in a [DISTANCE](DISTANCE.md) action before we computed the SECONDARY_STRUCTURE_RMSD values. These distances are then transformed by a Heavyside function that is 1 if
114 : the distance is less than 1 nm and zero otherwise. This vector of ones and zeros is then fed into the SECONDARY_STRUCTURE_RMSD action as a mask. By doing this we ensure that SECONDARY_STRUCTURE_RMSD
115 : values are only computed for parts of the strand that are close together. This avoids performing unnecessarily expensive calculations and is also a reasonable thing to do as if the two strands are more than
116 : 1 nm apart the residues in question are not at all similar to a beta sheet structure.
117 :
118 : */
119 : //+ENDPLUMEDOC
120 :
121 : namespace PLMD {
122 : namespace secondarystructure {
123 :
124 40 : class SecondaryStructureRMSDInput {
125 : /// The list of reference configurations
126 : std::vector<RMSD> myrmsd;
127 : public:
128 73 : static void registerKeywords( Keywords& keys ) {
129 73 : keys.setDisplayName("SECONDARY_STRUCTURE_RMSD");
130 73 : }
131 : static constexpr bool needsType=true;
132 : /// The number of atoms
133 : std::size_t natoms{0};
134 : /// The number of structures
135 : std::size_t nstructures{0};
136 : /// The number of indices per task
137 : std::size_t nindices_per_task{0};
138 : /// Are we operating without periodic boundary conditions
139 : bool nopbc{false};
140 : /// Variables for strands cutoff
141 : bool align_strands{false};
142 : /// The atoms involved in each of the secondary structure segments
143 : Matrix<unsigned> colvar_atoms;
144 : // void toACCDevice()const {
145 : // #pragma acc enter data copyin(this[0:1], myrmsd[0:nstructures],natoms,nstructures,nopbc,align_strands)
146 : // colvar_atoms.toACCDevice();
147 : // }
148 : // void removeFromACCDevice() const {
149 : // colvar_atoms.removeFromACCDevice();
150 : // #pragma acc exit data delete(align_strands,nopbc,nstructures,natoms,myrmsd[0:nstructures],this[0:1])
151 :
152 : static void calculateDistance( unsigned n, bool noderiv, const SecondaryStructureRMSDInput& actiondata, View<Vector> pos, ColvarOutput<double>& output );
153 : void setReferenceStructure( const std::string& type, double bondlength, std::vector<Vector>& structure );
154 80 : SecondaryStructureRMSDInput() = default;
155 : SecondaryStructureRMSDInput( const SecondaryStructureRMSDInput& m ) {
156 : *this = m;
157 : }
158 32 : SecondaryStructureRMSDInput& operator=( const SecondaryStructureRMSDInput& m ) {
159 32 : if( this==&m ) {
160 : return *this;
161 : }
162 32 : natoms = m.natoms;
163 32 : nstructures = m.nstructures;
164 32 : nindices_per_task = m.nindices_per_task;
165 32 : nopbc = m.nopbc;
166 32 : align_strands = m.align_strands;
167 32 : colvar_atoms=m.colvar_atoms;
168 : // setReferenceStructure appends to myrmsd, so the structures that are already here have to go
169 32 : myrmsd.clear();
170 79 : for(unsigned i=0; i<m.myrmsd.size(); ++i) {
171 47 : std::vector<Vector> mystruct( m.myrmsd[i].getReference() );
172 47 : setReferenceStructure( m.myrmsd[i].getMethod(), 0.0, mystruct );
173 47 : }
174 : return *this;
175 : }
176 : };
177 :
178 : typedef SecondaryStructureBase<SecondaryStructureRMSDInput> colv;
179 : PLUMED_REGISTER_ACTION(colv,"SECONDARY_STRUCTURE_RMSD")
180 :
181 96 : void SecondaryStructureRMSDInput::setReferenceStructure( const std::string& type, double bondlength, std::vector<Vector>& structure ) {
182 96 : if( !myrmsd.empty() && structure.size()!=natoms ) {
183 0 : plumed_merror("input structures have a different number of atoms");
184 : }
185 : // RMSD::set puts the center of the reference in the origin
186 96 : std::vector<double> align( structure.size(), 1.0 ), displace( structure.size(), 1.0 );
187 96 : RMSD newrmsd;
188 96 : newrmsd.clear();
189 96 : newrmsd.set(align,displace,structure,type,true,true);
190 96 : myrmsd.push_back( newrmsd );
191 96 : nstructures=myrmsd.size();
192 96 : natoms=structure.size();
193 192 : }
194 :
195 191632 : void SecondaryStructureRMSDInput::calculateDistance( unsigned n,
196 : bool noderiv,
197 : const SecondaryStructureRMSDInput& actiondata,
198 : const View<Vector> pos,
199 : ColvarOutput<double>& output ) {
200 191632 : std::vector<Vector> myderivs( actiondata.natoms );
201 191632 : output.values[n] = actiondata.myrmsd[n].calculate( pos, myderivs, false );
202 :
203 191632 : if( noderiv ) {
204 : return;
205 : }
206 92786 : Tensor v;
207 : v.zero();
208 2873522 : for(unsigned j=0; j<pos.size(); ++j) {
209 : output.derivs[n][j] = myderivs[j];
210 5561472 : v += (-1.0*Tensor( pos[j], myderivs[j] ));
211 : }
212 92786 : output.virial.set( n, v );
213 191632 : }
214 :
215 : }
216 : }
|
