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 4 : void AlphaBeta::registerKeywords(Keywords& keys) {
114 4 : ActionShortcut::registerKeywords( keys );
115 4 : keys.addFlag("NOPBC",false,"ignore the periodic boundary conditions when calculating distances");
116 4 : 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 8 : keys.reset_style("ATOMS","atoms");
120 4 : 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 4 : 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 8 : keys.setValueDescription("scalar","the alpha beta CV");
125 4 : keys.needsAction("CONSTANT");
126 4 : keys.needsAction("TORSION");
127 4 : keys.needsAction("COMBINE");
128 4 : keys.needsAction("CUSTOM");
129 4 : keys.needsAction("SUM");
130 4 : }
131 :
132 2 : AlphaBeta::AlphaBeta(const ActionOptions& ao):
133 : Action(ao),
134 2 : ActionShortcut(ao) {
135 : // Read in the reference value
136 : std::string refstr;
137 4 : parse("REFERENCE",refstr);
138 : unsigned nref=0;
139 2 : if( refstr.length()==0 ) {
140 : for(unsigned i=0;; ++i) {
141 : std::string refval;
142 18 : if( !parseNumbered( "REFERENCE", i+1, refval ) ) {
143 : break;
144 : }
145 8 : if( i==0 ) {
146 : refstr = refval;
147 : } else {
148 14 : refstr += "," + refval;
149 : }
150 8 : nref++;
151 8 : }
152 : }
153 : std::string coeffstr;
154 4 : parse("COEFFICIENT",coeffstr);
155 : unsigned ncoeff=0;
156 2 : if( coeffstr.length()==0 ) {
157 : for(unsigned i=0;; ++i) {
158 : std::string coeff;
159 4 : if( !parseNumbered( "COEFFICIENT", i+1, coeff) ) {
160 : break;
161 : }
162 0 : if( i==0 ) {
163 : coeffstr = coeff;
164 : } else {
165 0 : coeffstr += "," + coeff;
166 : }
167 0 : ncoeff++;
168 0 : }
169 : }
170 2 : if( coeffstr.length()==0 ) {
171 : coeffstr="1";
172 : }
173 : // Calculate angles
174 4 : readInputLine( getShortcutLabel() + "_torsions: TORSION " + convertInputLineToString() );
175 2 : ActionWithValue* av = plumed.getActionSet().selectWithLabel<ActionWithValue*>( getShortcutLabel() + "_torsions" );
176 2 : plumed_assert( av && (av->copyOutput(0))->getRank()==1 );
177 2 : if( nref==0 ) {
178 1 : std::string refval=refstr;
179 8 : for(unsigned i=1; i<(av->copyOutput(0))->getShape()[0]; ++i) {
180 14 : refstr += "," + refval;
181 : }
182 1 : } else if( nref!=(av->copyOutput(0))->getShape()[0] ) {
183 0 : error("mismatch between number of reference values and number of ATOMS specified");
184 : }
185 2 : if( ncoeff==0 ) {
186 2 : std::string coeff=coeffstr;
187 16 : for(unsigned i=1; i<(av->copyOutput(0))->getShape()[0]; ++i) {
188 28 : coeffstr += "," + coeff;
189 : }
190 0 : } else if( ncoeff!=(av->copyOutput(0))->getShape()[0] ) {
191 0 : error("mismatch between number of coefficients and number of ATOMS specified");
192 : }
193 4 : readInputLine( getShortcutLabel() + "_ref: CONSTANT VALUES=" + refstr );
194 4 : readInputLine( getShortcutLabel() + "_coeff: CONSTANT VALUES=" + coeffstr );
195 : // Caculate difference from reference using combine
196 4 : readInputLine( getShortcutLabel() + "_comb: COMBINE ARG=" + getShortcutLabel() + "_torsions," + getShortcutLabel() + "_ref COEFFICIENTS=1,-1 PERIODIC=NO" );
197 : // Now matheval for cosine bit
198 4 : readInputLine( getShortcutLabel() + "_cos: CUSTOM ARG=" + getShortcutLabel() + "_comb," + getShortcutLabel() + "_coeff FUNC=y*(0.5+0.5*cos(x)) PERIODIC=NO");
199 : // And combine to get final value
200 4 : readInputLine( getShortcutLabel() + ": SUM ARG=" + getShortcutLabel() + "_cos PERIODIC=NO");
201 2 : }
202 :
203 : }
204 : }
|
