Line data Source code
1 : /* +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
2 : Copyright (c) 2016-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 :
23 : /*
24 : This vast majority of the source code in this file was written by
25 : Sandro Bottaro with some help from Giovanni Bussi
26 : */
27 :
28 : #include "Colvar.h"
29 : #include "core/PlumedMain.h"
30 : #include "core/ActionRegister.h"
31 : #include "tools/PDB.h"
32 : #include "tools/ERMSD.h"
33 :
34 : namespace PLMD {
35 : namespace colvar {
36 :
37 :
38 : //+PLUMEDOC COLVAR ERMSD
39 : /*
40 : Calculate eRMSD with respect to a reference structure.
41 :
42 : eRMSD is a metric developed for measuring distances between three-dimensional RNA structures.
43 : The standard RMSD measure is highly inaccurate when measuring distances among three-dimensional
44 : structures of nucleic acids.
45 : It is not unusual, for example, that two RNA structures with low RMSD (i.e. less than 0.4nm) display a completely different network of base-base interactions.
46 :
47 : eRMSD measures the distance between structures by considering only the relative positions and orientations of nucleobases. The eRMSD can be considered as a vectorial version of contact maps and it is calculated as follows:
48 :
49 : 1. Set up a local reference system in the center of the six-membered ring of each nucleobase in a molecule.
50 : The xy plane lies on the plane of the nucleobase, and it is oriented such that the Watson-Crick interaction is always at $\theta\approx 60^{\circ}$.
51 :
52 : 2. Calculate all pairwise distance vectors $\vec{r}_{i,j}$ among base centers.
53 :
54 : 3. Rescale distance vectors as $\tilde{\vec{r}}_{i,j}=(r_x/a,r_y/a,r_z/b)$, where a=b=5 Angstrom, c=3 Angstrom. This rescaling has the effect of weighting more deviations on the z-axis with respect to the x/y directions.
55 :
56 : 4. Calculate the G vectors
57 :
58 : $$
59 : \vec{G}(\tilde{\vec{r}}) = (\sin(\gamma \tilde{r}) \tilde{r}_x/\tilde{r},\sin(\gamma \tilde{r}) \tilde{r}_y/\tilde{r},\sin(\gamma \tilde{r}) \tilde{r}_z/\tilde{r}, 1+\cos(\gamma \tilde{r})) \times
60 : \frac{\Theta(\tilde{r}_{cutoff}-\tilde{r})}{\gamma}
61 : $$
62 :
63 : Here, $ \gamma = \pi/\tilde{r}_{cutoff}$ and $ \Theta $ is the Heaviside step function. The default cutoff is set to 2.4.
64 :
65 : 5. The eRMSD between two structures $ \alpha $ and $ \beta $ reads
66 :
67 : $$
68 : eRMSD = \sqrt{\frac{1}{N} \sum_{j,k} \vert \vec{G}(\tilde{\vec{r}}_{jk}^{\alpha}) - \vec{G}(\tilde{\vec{r}}_{jk}^{\beta}) \vert^2 }
69 : $$
70 :
71 : Using the default cutoff, two structures with eRMSD of 0.7 or lower can be considered as significantly similar. A full description of the eRMSD can be found in the first paper in the bibliography below.
72 :
73 : ERMSD is computed using the position of three atoms on the 6-membered ring of each involved nucleobase. The atoms should be:
74 : - C2,C4,C6 for pyrimdines
75 : - C2,C6,C4 for purines
76 :
77 : The different order for purines and pyrimidines is fundamental and allows you to compute ERMSD between structures with different
78 : sequences as well! Notice that the simplest way to avoid mistakes in choosing these atoms is to use the `@lcs-#` strings
79 : as shown in the examples (see also [MOLINFO](MOLINFO.md)).
80 :
81 : !!! caution "ERMSD paths not possible"
82 :
83 : Notice that the ERMSD implemented here is not integrated with the other metrics in plumed. As a consequence, it is not (yet) possible
84 : to e.g. build path collective variables using ERMSD
85 :
86 : !!! caution "ERMSD with multiple molecules"
87 :
88 : Notice that ERMSD expect a single molecule and makes coordinates whole before anything else. As such, results might be unexpected
89 : for a multi molecular system.
90 :
91 : ## Examples
92 :
93 : Calculate the eRMSD from reference structure reference.pdb using the default cutoff (2.4). The list of residues involved in the calculation has to be specified. In this example, the eRMSD is calculated
94 : considering residues 1,2,3,4,5,6.
95 :
96 : ```plumed
97 : #SETTINGS MOLFILE=regtest/basic/rt-ermsd/ref.pdb
98 : MOLINFO STRUCTURE=regtest/basic/rt-ermsd/ref.pdb
99 : eRMSD1: ERMSD REFERENCE=regtest/basic/rt-ermsd/ref.pdb ATOMS=@lcs-1,@lcs-2,@lcs-3,@lcs-4,@lcs-5,@lcs-6
100 : ```
101 :
102 : */
103 : //+ENDPLUMEDOC
104 :
105 :
106 : class ERMSD : public Colvar {
107 :
108 :
109 : std::vector<Vector> derivs;
110 : PLMD::ERMSD ermsd;
111 : bool pbc;
112 :
113 : public:
114 : explicit ERMSD(const ActionOptions&);
115 : void calculate() override;
116 : static void registerKeywords(Keywords& keys);
117 : };
118 :
119 : PLUMED_REGISTER_ACTION(ERMSD,"ERMSD")
120 :
121 11 : void ERMSD::registerKeywords(Keywords& keys) {
122 11 : Colvar::registerKeywords(keys);
123 :
124 11 : keys.add("compulsory","REFERENCE","a file in pdb format containing the reference structure and the atoms involved in the CV.");
125 11 : keys.add("compulsory","CUTOFF","2.4","only pairs of atoms closer than CUTOFF are considered in the calculation.");
126 11 : keys.add("atoms","ATOMS","the list of atoms (use lcs)");
127 11 : keys.add("optional","PAIRS","List of pairs considered. All pairs are considered if this value is not specified.");
128 22 : keys.setValueDescription("scalar","the eRMSD between the instantaneous structure and the reference structure that was input");
129 11 : keys.addDOI("10.1093/nar/gku972");
130 11 : }
131 :
132 7 : ERMSD::ERMSD(const ActionOptions&ao):
133 7 : PLUMED_COLVAR_INIT(ao), pbc(true) {
134 : std::string reference;
135 7 : parse("REFERENCE",reference);
136 7 : double cutoff=2.4;
137 7 : parse("CUTOFF",cutoff);
138 :
139 :
140 7 : bool nopbc(false);
141 8 : parseFlag("NOPBC",nopbc);
142 7 : pbc=!nopbc;
143 :
144 7 : std::vector<AtomNumber> atoms_;
145 7 : parseAtomList("ATOMS",atoms_);
146 :
147 7 : std::vector<unsigned> pairs_;
148 7 : parseVector("PAIRS",pairs_);
149 7 : checkRead();
150 :
151 7 : addValueWithDerivatives();
152 7 : setNotPeriodic();
153 :
154 7 : if(atoms_.size()<6) {
155 0 : error("at least six atoms should be specified");
156 : }
157 7 : if(atoms_.size()%3!=0) {
158 0 : error("Atoms are not multiple of 3");
159 : }
160 7 : if(pairs_.size()%2!=0) {
161 0 : error("pairs are not multiple of 2");
162 : }
163 :
164 :
165 :
166 : // read everything in ang and transform to nm if we are not in natural units
167 7 : PDB pdb;
168 7 : if( !pdb.read(reference,usingNaturalUnits(),0.1/getUnits().getLength()) ) {
169 0 : error("missing input file " + reference );
170 : }
171 : // store target_ distance
172 7 : std::vector <Vector> reference_positions;
173 7 : unsigned natoms = atoms_.size();
174 7 : log.printf("Read %u atoms\n",natoms);
175 :
176 7 : reference_positions.resize(natoms);
177 877 : for(unsigned i=0; i<natoms; i++) {
178 870 : reference_positions[i] = pdb.getPosition(atoms_[i]);
179 : }
180 :
181 : // shift to count from zero (pairs are residue indices, 1-based in input)
182 32 : for(unsigned i=0; i<pairs_.size(); ++i) {
183 26 : if(pairs_[i]<1 || pairs_[i]>natoms/3) {
184 2 : error("PAIRS should contain residue indices between 1 and the number of residues (" + std::to_string(natoms/3) + ")");
185 : }
186 25 : pairs_[i]--;
187 : }
188 :
189 6 : ermsd.setReference(reference_positions,pairs_,cutoff/getUnits().getLength());
190 :
191 6 : requestAtoms(atoms_);
192 6 : derivs.resize(natoms);
193 :
194 6 : log.printf(" reference from file %s\n",reference.c_str());
195 6 : log.printf(" which contains %u atoms\n",natoms);
196 :
197 6 : log<<" Bibliography "
198 6 : <<plumed.cite("Bottaro, Di Palma, and Bussi, Nucleic Acids Res. 42, 13306 (2014)")
199 18 : <<plumed.cite("Bottaro, Banas, Sponer, and Bussi, J. Phys. Chem. Lett. 7, 4032 (2016)")<<"\n";
200 :
201 19 : }
202 :
203 : // calculator
204 228 : void ERMSD::calculate() {
205 : // set derivatives to zero
206 228 : Tools::set_to_zero(derivs);
207 : double ermsdist;
208 228 : Tensor virial;
209 : // This is a trick to avoid explicit virial calculation
210 : // 1. we make the molecule whole (unless NOPBC was requested)
211 228 : if(pbc) {
212 225 : makeWhole();
213 : }
214 : // 2. we ignore pbcs
215 228 : Pbc fake_pbc;
216 : // Notice that this might have problems when having 2 RNA molecules (hybridization).
217 :
218 228 : ermsdist=ermsd.calculate(getPositions(),fake_pbc,derivs,virial);
219 228 : const double scale=getUnits().getLength();
220 228 : setValue(ermsdist*scale);
221 :
222 47550 : for(unsigned i=0; i<derivs.size(); ++i) {
223 47322 : setAtomsDerivatives(i,derivs[i]*scale);
224 : }
225 :
226 228 : setBoxDerivativesNoPbc();
227 228 : }
228 :
229 : }
230 : }
|
