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 : #ifdef __PLUMED_HAS_OPENACC
23 : #define __PLUMED_USE_OPENACC 1
24 : #endif //__PLUMED_HAS_OPENACC
25 : #include "SecondaryStructureBase.h"
26 : #include "core/ActionRegister.h"
27 : #include "tools/Matrix.h"
28 :
29 : #include<vector>
30 :
31 : //+PLUMEDOC MCOLVAR SECONDARY_STRUCTURE_DRMSD
32 : /*
33 : Calclulate the DRMSD between segments of a protein and a reference structure of interest
34 :
35 : This action is used in the shortcuts [ALPHARMSD](ALPHARMSD.md), [ANTIBETARMSD](ANTIBETARMSD.md) and [PARABETARMSD](PARABETARMSD.md). It calculates a
36 : vector of [DRMSD](DRMSD.md) values between a single reference multiple configurations and the instantaneous
37 : positions of various groups of atoms. For example, in the following input we define a single set of reference
38 : set of coordinates for 3 atoms.
39 :
40 : ```plumed
41 : c1: SECONDARY_STRUCTURE_DRMSD ...
42 : ATOMS=13-24
43 : BONDLENGTH=0.17 STRUCTURE1=1,0,0,0,1,0,0,0,1
44 : SEGMENT1=0,1,2 SEGMENT2=3,4,5
45 : SEGMENT3=6,7,8 SEGMENT4=9,10,11
46 : ...
47 : PRINT ARG=c1 FILE=colvar
48 : ```
49 :
50 : A four dimensional vector is then returned that contains the DRMSD distances between the 4 sets of 3 atoms that were specified using the `SEGMENT` keywords
51 : and the reference coordinates. Notice that you can use multiple instances of the `STRUCTURE` keyword. In general the the number of vectors output
52 : 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
53 : of atoms specified by the ATOMS keyword. `SEGMENT1=0,1,2` in the above input is measuring the distance between the instaneous positions of atoms 13, 14 and 15
54 : and the reference structure.
55 :
56 : ## Periodic boundary conditions
57 :
58 : If you turn off periodic boundary conditions by using the NOPBC flag as shown below:
59 :
60 : ```plumed
61 : c1: SECONDARY_STRUCTURE_DRMSD ...
62 : ATOMS=13-24 NOPBC
63 : BONDLENGTH=0.17 STRUCTURE1=1,0,0,0,1,0,0,0,1
64 : SEGMENT1=0,1,2 SEGMENT2=3,4,5
65 : SEGMENT3=6,7,8 SEGMENT4=9,10,11
66 : ...
67 : PRINT ARG=c1 FILE=colvar
68 : ```
69 :
70 : the distances between in the instaneous structure are evaluated in a way that does not take the periodic boundary conditions
71 : into account. Whenver this flag is __not__ present the distances in the instaneous structure are evaluated in a way that takes the periodic boundary conditions
72 : into account.
73 :
74 : If you use the `ALIGN_STRANDS` flag to evaluate the distances in a way that takes the periodic boundary conditions into accounts means that
75 : calculating the distance between the positions of the 7th and 22nd atoms ($\mathbf{x}_7$ and $\mathbf{x}_{22}$) of the instanenous structures
76 : using the periodic boundary conditions. The 22 atom is the repositioned at:
77 :
78 : $$
79 : \mathbf{x}_{22} = \mathbf{x}_7 + ( \mathbf{x}_{22} - \mathbf{x}_7 )
80 : $$
81 :
82 : 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
83 : 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 perfomed and the
84 : 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)
85 : sheets are not broken by the periodic boundaries.
86 :
87 : ## Working with beta sheet like structures
88 :
89 : As discussed in the documentation for [ANTIBETARMSD](ANTIBETARMSD.md) and [PARABETARMSD](PARABETARMSD.md) using the STRANDS_KEYWORD cutoff can speed up the calculation
90 : dramatically. The reason this keyword gives such dramatric speed ups is illustrated by the example input:
91 :
92 : ```plumed
93 : ab_cut_dists: DISTANCE ...
94 : ATOMS1=19,69 ATOMS2=19,79
95 : ATOMS3=19,89 ATOMS4=19,99
96 : ATOMS5=19,109
97 : ...
98 : ab_cut: CUSTOM ARG=ab_cut_dists FUNC=step(1-x) PERIODIC=NO
99 : ab_rmsd: SECONDARY_STRUCTURE_DRMSD ...
100 : BONDLENGTH=0.17 ALIGN_STRANDS MASK=ab_cut
101 : 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
102 : 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
103 : 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
104 : 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
105 : 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
106 : 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
107 : 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
108 : ...
109 : ab_ltu: LESS_THAN ARG=ab_rmsd SWITCH={RATIONAL R_0=0.1 D_0=0.0 NN=8 MM=12} MASK=ab_cut
110 : ab_lt: CUSTOM ARG=ab_ltu,ab_cut FUNC=x*y PERIODIC=NO
111 : ab: SUM ARG=ab_lt PERIODIC=NO
112 : PRINT ARG=ab FILE=colvar
113 : ```
114 :
115 : 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
116 : beta sheet are computed in a [DISTANCE](DISTANCE.md) action before we computed the SECONDARY_STRUCTURE_DRMSD values. These distances are then transformed by a Heavyside function that is 1 if
117 : the distance is less than 1 nm and zero otherwise. This vector of ones and zeros is then fed into the SECONDARY_STRUCTURE_DRMSD action as a mask. By doing this we ensure that SECONDARY_STRUCTURE_DRMSD
118 : values are only computed for parts of the strand that are close together. This avoids performing unecessarily expensive calculations and is also a reasonable thing to do as if the two strands are more than
119 : 1 nm appart the residues in question are not at all similar to a beta sheet structure.
120 :
121 : */
122 : //+ENDPLUMEDOC
123 :
124 : namespace PLMD {
125 : namespace secondarystructure {
126 :
127 : ///This is a sort of compact std::vector<std::vector<T>>
128 : template<typename T>
129 : class flexibleMemory {
130 : ///Stores the data
131 : std::vector<T> data{};
132 : ///Stores the boundaries
133 : std::vector<size_t> sizes{};
134 : //helpers for openACC
135 : T* d;
136 : size_t* sz;
137 : void update() {
138 84 : d=data.data();
139 84 : sz=sizes.data();
140 34 : }
141 : public:
142 : flexibleMemory() {
143 : update();
144 : };
145 : ~flexibleMemory() =default;
146 : flexibleMemory(const flexibleMemory& m) :
147 : data(m.data),
148 : sizes(m.sizes) {
149 : update();
150 : }
151 34 : flexibleMemory& operator=(const flexibleMemory& m) {
152 34 : if (this!=&m) {
153 34 : data=m.data;
154 34 : sizes=m.sizes;
155 : update();
156 : }
157 34 : return *this;
158 : }
159 : constexpr size_t size() const {
160 : return sizes.size();
161 : }
162 : //Access the i-th view of the data, use this to access the elements
163 169823 : constexpr View< const T> get(size_t i) const {
164 : //use this at the start of your function, to not calculate again and again the initial index
165 : size_t init=0;
166 231298 : for(size_t j=0; j<i; j++) {
167 61475 : init+=sz[j];
168 : }
169 169823 : return View<const T>(d+init,sz[i]);
170 : }
171 50 : void extend(const std::vector<T>& v) {
172 50 : data.insert(data.end(),v.begin(),v.end());
173 50 : sizes.push_back(v.size());
174 : update();
175 50 : }
176 : void toACCDevice()const {
177 : #pragma acc enter data copyin(this[0:1], d[0:data.size()], sz[0:sizes.size()])
178 : }
179 : void removeFromACCDevice() const {
180 : #pragma acc exit data delete(sz[0:sizes.size()], d[0:data.size()], this[0:1])
181 : }
182 : };
183 :
184 : class SecondaryStructureDRMSDInput {
185 : private:
186 : struct pairDistance {
187 : double distance;
188 : unsigned i;
189 : unsigned j;
190 : };
191 : /// The list of reference configurations
192 : flexibleMemory<pairDistance> drmsd_targets{};
193 : flexibleMemory<unsigned> drmsd_atoms{};
194 : /// The general input for the secondary structure variable
195 : public:
196 : size_t natoms{0};
197 : size_t nstructures{0};
198 : size_t nindices_per_task{0};
199 : /// Are we operating without periodic boundary conditions
200 : bool nopbc{false};
201 : /// Variables for strands cutoff
202 : bool align_strands{false};
203 : /// The atoms involved in each of the secondary structure segments
204 : Matrix<unsigned> colvar_atoms{};
205 : static void calculateDistance( unsigned n,
206 : bool noderiv,
207 : const SecondaryStructureDRMSDInput& actiondata,
208 : View<Vector> pos,
209 : ColvarOutput& output );
210 : void setReferenceStructure( std::string type, double bondlength, std::vector<Vector>& structure );
211 : void toACCDevice()const {
212 : #pragma acc enter data copyin(this[0:1], natoms,nstructures,nindices_per_task,nopbc,align_strands)
213 : drmsd_targets.toACCDevice();
214 : drmsd_atoms.toACCDevice();
215 : colvar_atoms.toACCDevice();
216 : }
217 : void removeFromACCDevice() const {
218 : colvar_atoms.removeFromACCDevice();
219 : drmsd_atoms.removeFromACCDevice();
220 : drmsd_targets.removeFromACCDevice();
221 : #pragma acc exit data delete(align_strands,nopbc,nindices_per_task,nstructures,natoms,this[0:1])
222 :
223 : }
224 : };
225 :
226 : typedef SecondaryStructureBase<SecondaryStructureDRMSDInput> colv;
227 : PLUMED_REGISTER_ACTION(colv,"SECONDARY_STRUCTURE_DRMSD")
228 :
229 25 : void SecondaryStructureDRMSDInput::setReferenceStructure(
230 : std::string /*type*/,
231 : double bondlength,
232 : std::vector<Vector>& structure ) {
233 : std::vector<pairDistance> targets;
234 : //set is ordered and contains no duplicated data
235 : std::set<unsigned> atoms_targets;
236 : //targets.reserve( (structure.size()*(structure.size()-1))/2 );?
237 750 : for(unsigned i=0; i<structure.size()-1; ++i) {
238 11600 : for(unsigned j=i+1; j<structure.size(); ++j) {
239 10875 : double distance = delta( structure[i], structure[j] ).modulo();
240 10875 : if(distance > bondlength) {
241 10172 : targets.emplace_back(pairDistance{distance,i,j});
242 10172 : atoms_targets.emplace(i);
243 10172 : atoms_targets.emplace(j);
244 : }
245 : }
246 : }
247 25 : drmsd_targets.extend( targets );
248 :
249 50 : drmsd_atoms.extend(std::vector<unsigned> {atoms_targets.begin(),atoms_targets.end()});
250 25 : nstructures = drmsd_targets.size();
251 25 : if (natoms==0) {
252 17 : natoms = structure.size();
253 8 : } else if (natoms!=structure.size()) {
254 0 : plumed_merror("input structures have a different number of atoms");
255 : }
256 :
257 25 : }
258 :
259 68042 : void SecondaryStructureDRMSDInput::calculateDistance(
260 : const unsigned n,
261 : const bool noderiv,
262 : const SecondaryStructureDRMSDInput& actiondata,
263 : const View<Vector> pos,
264 : ColvarOutput& output ) {
265 : {
266 68042 : const auto targetList = actiondata.drmsd_atoms.get(n);
267 : const auto targetAtoms=targetList.size();
268 68042 : if( !noderiv ) {
269 33739 : output.virial.set( n, Tensor(0,0,0,0,0,0,0,0,0) );
270 1045909 : for(unsigned i=0; i<targetAtoms; ++i ) {
271 1012170 : output.derivs[n][targetList[i]] =Vector(0.0,0.0,0.0);
272 : }
273 : }
274 : }
275 :
276 : double drmsd = 0.0;
277 : Vector distance;
278 : Vector dd;
279 :
280 68042 : const auto target = actiondata.drmsd_targets.get(n);
281 : const auto targetSize=target.size();
282 :
283 27760860 : for (unsigned i=0; i<targetSize; ++i) {
284 82848990 : const auto& [reference,k,j] = target[i];
285 27692818 : distance=delta( pos[k], pos[j] );
286 27692818 : const double len = distance.modulo();
287 27692818 : const double diff = len - reference;
288 27692818 : drmsd += diff*diff;
289 :
290 27692818 : if( !noderiv ) {
291 13731677 : dd=distance*(diff / len);
292 13731677 : output.derivs[n][k] -= dd;//-3 mul
293 13731677 : output.derivs[n][j] += dd;//-3mul
294 13731677 : output.virial.set( n, output.virial[n] - Tensor(dd,distance) );//-9 mul
295 : }
296 : }
297 68042 : const double inpairs = 1./static_cast<double>(targetSize);
298 68042 : output.values[n] = sqrt(inpairs*drmsd);
299 :
300 68042 : if( !noderiv ) {
301 33739 : double scalef = inpairs / output.values[n];
302 :
303 33739 : PLMD::LoopUnroller<9>::_mul(output.virial.getView(n).data(),scalef);
304 33739 : const auto targetList = actiondata.drmsd_atoms.get(n);
305 1045909 : for(unsigned i=0; i<actiondata.natoms; ++i ) {
306 1012170 : output.derivs[n][targetList[i]] *= scalef;
307 : }
308 : }
309 68042 : }
310 :
311 : }
312 : }
|
