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 "core/ActionShortcut.h"
23 : #include "core/ActionWithValue.h"
24 : #include "core/ActionRegister.h"
25 : #include "core/PlumedMain.h"
26 : #include "core/ActionSet.h"
27 : #include <string>
28 : #include <cmath>
29 :
30 : //+PLUMEDOC MCOLVAR ALPHABETA
31 : /*
32 : Measures a distance including pbc between the instantaneous values of a set of torsional angles and set of reference values.
33 :
34 : This shortcut calculates the following quantity.
35 :
36 : $$
37 : s = \frac{1}{2} \sum_i c_i \left[ 1 + \cos( \phi_i - \phi_i^{\textrm{Ref}} ) \right]
38 : $$
39 :
40 : where the $\phi_i$ values are the instantaneous values for the [TORSION](TORSION.md) angles of interest.
41 : The $\phi_i^{\textrm{Ref}}$ values are reference values for the torsional angles that are specified in the input file.
42 :
43 : The following provides an example of the input for an alpha beta similarity.
44 :
45 : ```plumed
46 : ab: ALPHABETA ...
47 : ATOMS1=168,170,172,188 REFERENCE1=3.14
48 : ATOMS2=170,172,188,190 REFERENCE2=3.14
49 : ATOMS3=188,190,192,230 REFERENCE3=3.14
50 : ...
51 : PRINT ARG=ab FILE=colvar STRIDE=10
52 : ```
53 :
54 : Because all the reference values are the same we can also calculate the same quantity using
55 :
56 : ```plumed
57 : ab: ALPHABETA ...
58 : ATOMS1=168,170,172,188 REFERENCE=3.14
59 : ATOMS2=170,172,188,190
60 : ATOMS3=188,190,192,230
61 : ...
62 : PRINT ARG=ab FILE=colvar STRIDE=10
63 : ```
64 :
65 : Writing out the atoms involved in all the torsion angles in this way can be rather tedious. Thankfully if you are working with protein you
66 : can avoid this by using the [MOLINFO](MOLINFO.md) command. PLUMED uses the pdb file that you provide to this command to learn
67 : about the topology of the protein molecule. This means that you can specify torsion angles using the following syntax:
68 :
69 : ```plumed
70 : #SETTINGS MOLFILE=regtest/basic/rt32/helix.pdb
71 : MOLINFO MOLTYPE=protein STRUCTURE=regtest/basic/rt32/helix.pdb
72 : ab: ALPHABETA ...
73 : ATOMS1=@phi-3 REFERENCE=3.14 COEFFICIENT1=2
74 : ATOMS2=@psi-3 COEFFICIENT2=0.5
75 : ATOMS3=@phi-4 COEFFICIENT3=1
76 : ...
77 : PRINT ARG=ab FILE=colvar STRIDE=10
78 : ```
79 :
80 : Here, `@phi-3` tells plumed that you would like to calculate the $\phi$ angle in the third residue of the protein.
81 : Similarly `@psi-4` tells plumed that you want to calculate the $\psi$ angle of the fourth residue of the protein.
82 : Notice, also, that in the first two examples the coefficients $c_i$ in the expression above were all set equal to one.
83 : In the example above we use the COEFFICIENT keywords to set these quantities to three different values.
84 :
85 : Notice, last of all, that in the above examples we reassemble any molecules that have been broken by the periodic boundary
86 : conditions using a procedure like that used in [WHOLEMOLECULES](WHOLEMOLECULES.md) before calculating the torsion angles.
87 : If you wish to turn this off for any reason you use the NOPBC flag as shown below:
88 :
89 : ```plumed
90 : ab: ALPHABETA ...
91 : ATOMS1=168,170,172,188 REFERENCE=3.14
92 : ATOMS2=170,172,188,190 NOPBC
93 : ATOMS3=188,190,192,230
94 : ...
95 : PRINT ARG=ab FILE=colvar STRIDE=10
96 : ```
97 :
98 :
99 : */
100 : //+ENDPLUMEDOC
101 :
102 : namespace PLMD {
103 : namespace multicolvar {
104 :
105 : class AlphaBeta : public ActionShortcut {
106 : public:
107 : static void registerKeywords(Keywords& keys);
108 : explicit AlphaBeta(const ActionOptions&);
109 : };
110 :
111 : PLUMED_REGISTER_ACTION(AlphaBeta,"ALPHABETA")
112 :
113 7 : void AlphaBeta::registerKeywords(Keywords& keys) {
114 7 : ActionShortcut::registerKeywords( keys );
115 7 : keys.addFlag("NOPBC",false,"ignore the periodic boundary conditions when calculating distances");
116 7 : keys.add("numbered","ATOMS","the atoms involved for each of the torsions you wish to calculate. "
117 : "Keywords like ATOMS1, ATOMS2, ATOMS3,... should be listed and one torsion will be "
118 : "calculated for each ATOM keyword you specify");
119 14 : keys.reset_style("ATOMS","atoms");
120 7 : keys.add("numbered","REFERENCE","the reference values for each of the torsional angles. If you use a single REFERENCE value the "
121 : "same reference value is used for all torsions");
122 7 : keys.add("numbered","COEFFICIENT","the coefficient for each of the torsional angles. If you use a single COEFFICIENT value the "
123 : "same reference value is used for all torsional angles");
124 14 : keys.setValueDescription("scalar","the alpha beta CV");
125 7 : keys.needsAction("CONSTANT");
126 7 : keys.needsAction("TORSION");
127 7 : keys.needsAction("COMBINE");
128 7 : keys.needsAction("CUSTOM");
129 7 : keys.needsAction("SUM");
130 7 : }
131 :
132 3 : AlphaBeta::AlphaBeta(const ActionOptions& ao):
133 : Action(ao),
134 3 : ActionShortcut(ao) {
135 : // Read in the reference value
136 : std::string refstr;
137 6 : parse("REFERENCE",refstr);
138 : unsigned nref=0;
139 3 : if( refstr.length()==0 ) {
140 : for(unsigned i=0;; ++i) {
141 : std::string refval;
142 20 : if( !parseNumbered( "REFERENCE", i+1, refval ) ) {
143 : break;
144 : }
145 8 : if( i==0 ) {
146 : refstr = refval;
147 : } else {
148 7 : refstr += "," + refval;
149 : }
150 8 : nref++;
151 8 : }
152 : }
153 3 : if( refstr.length()==0 ) {
154 2 : error("no reference values found: use REFERENCE to give the same value for all the torsions or REFERENCE1, REFERENCE2, ... to give one for each");
155 : }
156 : std::string coeffstr;
157 4 : parse("COEFFICIENT",coeffstr);
158 : unsigned ncoeff=0;
159 2 : if( coeffstr.length()==0 ) {
160 : for(unsigned i=0;; ++i) {
161 : std::string coeff;
162 4 : if( !parseNumbered( "COEFFICIENT", i+1, coeff) ) {
163 : break;
164 : }
165 0 : if( i==0 ) {
166 : coeffstr = coeff;
167 : } else {
168 0 : coeffstr += "," + coeff;
169 : }
170 0 : ncoeff++;
171 0 : }
172 : }
173 2 : if( coeffstr.length()==0 ) {
174 : coeffstr="1";
175 : }
176 : // Calculate angles
177 4 : readInputLine( getShortcutLabel() + "_torsions: TORSION " + convertInputLineToString() );
178 2 : ActionWithValue* av = plumed.getActionSet().selectWithLabel<ActionWithValue*>( getShortcutLabel() + "_torsions" );
179 2 : plumed_assert( av && (av->copyOutput(0))->getRank()==1 );
180 2 : if( nref==0 ) {
181 1 : std::string refval=refstr;
182 8 : for(unsigned i=1; i<(av->copyOutput(0))->getShape()[0]; ++i) {
183 7 : refstr += "," + refval;
184 : }
185 1 : } else if( nref!=(av->copyOutput(0))->getShape()[0] ) {
186 0 : error("mismatch between number of reference values and number of ATOMS specified");
187 : }
188 2 : if( ncoeff==0 ) {
189 2 : std::string coeff=coeffstr;
190 16 : for(unsigned i=1; i<(av->copyOutput(0))->getShape()[0]; ++i) {
191 14 : coeffstr += "," + coeff;
192 : }
193 0 : } else if( ncoeff!=(av->copyOutput(0))->getShape()[0] ) {
194 0 : error("mismatch between number of coefficients and number of ATOMS specified");
195 : }
196 2 : readInputLine( getShortcutLabel() + "_ref: CONSTANT VALUES=" + refstr );
197 2 : readInputLine( getShortcutLabel() + "_coeff: CONSTANT VALUES=" + coeffstr );
198 : // Caculate difference from reference using combine
199 6 : readInputLine( getShortcutLabel() + "_comb: COMBINE ARG=" + getShortcutLabel() + "_torsions," + getShortcutLabel() + "_ref COEFFICIENTS=1,-1 PERIODIC=NO" );
200 : // Now matheval for cosine bit
201 6 : readInputLine( getShortcutLabel() + "_cos: CUSTOM ARG=" + getShortcutLabel() + "_comb," + getShortcutLabel() + "_coeff FUNC=y*(0.5+0.5*cos(x)) PERIODIC=NO");
202 : // And combine to get final value
203 6 : readInputLine( getShortcutLabel() + ": SUM ARG=" + getShortcutLabel() + "_cos PERIODIC=NO");
204 32 : }
205 :
206 : }
207 : }
|
