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/RMSD.h"
28 :
29 : //+PLUMEDOC MCOLVAR SECONDARY_STRUCTURE_RMSD
30 : /*
31 : Calclulate the RMSD between segments of a protein and a reference structure of interest
32 :
33 : This action is used in the shortcuts [ALPHARMSD](ALPHARMSD.md), [ANTIBETARMSD](ANTIBETARMSD.md) and [PARABETARMSD](PARABETARMSD.md). It calculates a
34 : vector of [RMSD](RMSD.md) values between a single reference multiple configurations and the instantaneous
35 : positions of various groups of atoms. For example, in the following input we define a single set of reference
36 : set of coordinates for 3 atoms.
37 :
38 : ```plumed
39 : c1: SECONDARY_STRUCTURE_RMSD ...
40 : ATOMS=13-24
41 : STRUCTURE1=1,0,0,0,1,0,0,0,1
42 : SEGMENT1=0,1,2 SEGMENT2=3,4,5
43 : SEGMENT3=6,7,8 SEGMENT4=9,10,11
44 : TYPE=OPTIMAL
45 : ...
46 : PRINT ARG=c1 FILE=colvar
47 : ```
48 :
49 : 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
50 : and the reference coordinates. Notice that you can use multiple instances of the `STRUCTURE` keyword. In general the the number of vectors output
51 : 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
52 : 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
53 : and the reference structure.
54 :
55 : ## Periodic boundary conditions
56 :
57 : If you turn off periodic boundary conditions by using the NOPBC flag as shown below:
58 :
59 : ```plumed
60 : c1: SECONDARY_STRUCTURE_RMSD ...
61 : ATOMS=13-24 NOPBC
62 : STRUCTURE1=1,0,0,0,1,0,0,0,1
63 : SEGMENT1=0,1,2 SEGMENT2=3,4,5
64 : SEGMENT3=6,7,8 SEGMENT4=9,10,11
65 : TYPE=SIMPLE
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_RMSD ...
100 : TYPE=OPTIMAL 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_RMSD 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_RMSD action as a mask. By doing this we ensure that SECONDARY_STRUCTURE_RMSD
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 21 : class SecondaryStructureRMSDInput {
128 : /// The list of reference configurations
129 : std::vector<RMSD> myrmsd;
130 : public:
131 : /// The number of atoms
132 : std::size_t natoms;
133 : /// The number of structures
134 : std::size_t nstructures;
135 : /// The number of indices per task
136 : std::size_t nindices_per_task;
137 : /// Are we operating without periodic boundary conditions
138 : bool nopbc;
139 : /// Variables for strands cutoff
140 : bool align_strands;
141 : /// The atoms involved in each of the secondary structure segments
142 : Matrix<unsigned> colvar_atoms;
143 : // void toACCDevice()const {
144 : // #pragma acc enter data copyin(this[0:1], myrmsd[0:nstructures],natoms,nstructures,nopbc,align_strands)
145 : // colvar_atoms.toACCDevice();
146 : // }
147 : // void removeFromACCDevice() const {
148 : // colvar_atoms.removeFromACCDevice();
149 : // #pragma acc exit data delete(align_strands,nopbc,nstructures,natoms,myrmsd[0:nstructures],this[0:1])
150 :
151 : // }
152 : static void calculateDistance( unsigned n, bool noderiv, const SecondaryStructureRMSDInput& actiondata, View<Vector> pos, ColvarOutput& output );
153 : void setReferenceStructure( const std::string& type, double bondlength, std::vector<Vector>& structure );
154 21 : SecondaryStructureRMSDInput& operator=( const SecondaryStructureRMSDInput& m ) {
155 21 : natoms = m.natoms;
156 21 : nstructures = m.nstructures;
157 21 : nindices_per_task = m.nindices_per_task;
158 21 : nopbc = m.nopbc;
159 21 : align_strands = m.align_strands;
160 21 : colvar_atoms=m.colvar_atoms;
161 52 : for(unsigned i=0; i<m.myrmsd.size(); ++i) {
162 31 : std::vector<Vector> mystruct( m.myrmsd[i].getReference() );
163 62 : setReferenceStructure( m.myrmsd[i].getMethod(), 0.0, mystruct );
164 : }
165 21 : return *this;
166 : }
167 : };
168 :
169 : typedef SecondaryStructureBase<SecondaryStructureRMSDInput> colv;
170 : PLUMED_REGISTER_ACTION(colv,"SECONDARY_STRUCTURE_RMSD")
171 :
172 62 : void SecondaryStructureRMSDInput::setReferenceStructure( const std::string& type, double bondlength, std::vector<Vector>& structure ) {
173 : Vector center;
174 62 : std::vector<double> align( structure.size(), 1.0 ), displace( structure.size(), 1.0 );
175 1922 : for(unsigned i=0; i<structure.size(); ++i) {
176 1860 : center+=structure[i]*align[i];
177 : }
178 1922 : for(unsigned i=0; i<structure.size(); ++i) {
179 : structure[i] -= center;
180 : }
181 62 : RMSD newrmsd;
182 62 : newrmsd.clear();
183 62 : newrmsd.set(align,displace,structure,type,true,true);
184 62 : myrmsd.push_back( newrmsd );
185 62 : nstructures=myrmsd.size();
186 62 : natoms=structure.size();
187 124 : }
188 :
189 145592 : void SecondaryStructureRMSDInput::calculateDistance( unsigned n,
190 : bool noderiv,
191 : const SecondaryStructureRMSDInput& actiondata,
192 : const View<Vector> pos,
193 : ColvarOutput& output ) {
194 145592 : std::vector<Vector> myderivs( actiondata.natoms );
195 145592 : output.values[n] = actiondata.myrmsd[n].calculate( pos, myderivs, false );
196 :
197 145592 : if( noderiv ) {
198 : return;
199 : }
200 : Tensor v;
201 : v.zero();
202 2256676 : for(unsigned j=0; j<pos.size(); ++j) {
203 : output.derivs[n][j] = myderivs[j];
204 4367760 : v += (-1.0*Tensor( pos[j], myderivs[j] ));
205 : }
206 72796 : output.virial.set( n, v );
207 : }
208 :
209 : }
210 : }
|
