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 \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 : 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 : LABEL=ab
51 : ... ALPHABETA
52 : PRINT ARG=ab FILE=colvar STRIDE=10
53 : ```
54 :
55 : Because all the reference values are the same we can also calculate the same quantity using
56 :
57 : ```plumed
58 : ALPHABETA ...
59 : ATOMS1=168,170,172,188 REFERENCE=3.14
60 : ATOMS2=170,172,188,190
61 : ATOMS3=188,190,192,230
62 : LABEL=ab
63 : ... ALPHABETA
64 : PRINT ARG=ab FILE=colvar STRIDE=10
65 : ```
66 :
67 : 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
68 : can avoid this by using the [MOLINFO](MOLINFO.md) command. PLUMED uses the pdb file that you provide to this command to learn
69 : about the topology of the protein molecule. This means that you can specify torsion angles using the following syntax:
70 :
71 : ```plumed
72 : #SETTINGS MOLFILE=regtest/basic/rt32/helix.pdb
73 : MOLINFO MOLTYPE=protein STRUCTURE=regtest/basic/rt32/helix.pdb
74 : ALPHABETA ...
75 : ATOMS1=@phi-3 REFERENCE=3.14
76 : ATOMS2=@psi-3
77 : ATOMS3=@phi-4
78 : LABEL=ab
79 : ... ALPHABETA
80 : PRINT ARG=ab FILE=colvar STRIDE=10
81 : ```
82 :
83 : Here, `@phi-3` tells plumed that you would like to calculate the $\phi$ angle in the third residue of the protein.
84 : Similarly `@psi-4` tells plumed that you want to calculate the $\psi$ angle of the fourth residue of the protein.
85 :
86 : */
87 : //+ENDPLUMEDOC
88 :
89 : namespace PLMD {
90 : namespace multicolvar {
91 :
92 : class AlphaBeta : public ActionShortcut {
93 : public:
94 : static void registerKeywords(Keywords& keys);
95 : explicit AlphaBeta(const ActionOptions&);
96 : };
97 :
98 : PLUMED_REGISTER_ACTION(AlphaBeta,"ALPHABETA")
99 :
100 4 : void AlphaBeta::registerKeywords(Keywords& keys) {
101 4 : ActionShortcut::registerKeywords( keys );
102 4 : keys.addFlag("NOPBC",false,"ignore the periodic boundary conditions when calculating distances");
103 4 : keys.add("numbered","ATOMS","the atoms involved for each of the torsions you wish to calculate. "
104 : "Keywords like ATOMS1, ATOMS2, ATOMS3,... should be listed and one torsion will be "
105 : "calculated for each ATOM keyword you specify");
106 8 : keys.reset_style("ATOMS","atoms");
107 4 : keys.add("numbered","REFERENCE","the reference values for each of the torsional angles. If you use a single REFERENCE value the "
108 : "same reference value is used for all torsions");
109 4 : keys.add("numbered","COEFFICIENT","the coefficient for each of the torsional angles. If you use a single COEFFICIENT value the "
110 : "same reference value is used for all torsional angles");
111 8 : keys.setValueDescription("scalar","the alpha beta CV");
112 4 : keys.needsAction("CONSTANT");
113 4 : keys.needsAction("TORSION");
114 4 : keys.needsAction("COMBINE");
115 4 : keys.needsAction("CUSTOM");
116 4 : keys.needsAction("SUM");
117 4 : }
118 :
119 2 : AlphaBeta::AlphaBeta(const ActionOptions& ao):
120 : Action(ao),
121 2 : ActionShortcut(ao) {
122 : // Read in the reference value
123 : std::string refstr;
124 4 : parse("REFERENCE",refstr);
125 : unsigned nref=0;
126 2 : if( refstr.length()==0 ) {
127 : for(unsigned i=0;; ++i) {
128 : std::string refval;
129 18 : if( !parseNumbered( "REFERENCE", i+1, refval ) ) {
130 : break;
131 : }
132 8 : if( i==0 ) {
133 : refstr = refval;
134 : } else {
135 14 : refstr += "," + refval;
136 : }
137 8 : nref++;
138 8 : }
139 : }
140 : std::string coeffstr;
141 4 : parse("COEFFICIENT",coeffstr);
142 : unsigned ncoeff=0;
143 2 : if( coeffstr.length()==0 ) {
144 : for(unsigned i=0;; ++i) {
145 : std::string coeff;
146 4 : if( !parseNumbered( "COEFFICIENT", i+1, coeff) ) {
147 : break;
148 : }
149 0 : if( i==0 ) {
150 : coeffstr = coeff;
151 : } else {
152 0 : coeffstr += "," + coeff;
153 : }
154 0 : ncoeff++;
155 0 : }
156 : }
157 2 : if( coeffstr.length()==0 ) {
158 : coeffstr="1";
159 : }
160 : // Calculate angles
161 4 : readInputLine( getShortcutLabel() + "_torsions: TORSION " + convertInputLineToString() );
162 2 : ActionWithValue* av = plumed.getActionSet().selectWithLabel<ActionWithValue*>( getShortcutLabel() + "_torsions" );
163 2 : plumed_assert( av && (av->copyOutput(0))->getRank()==1 );
164 2 : if( nref==0 ) {
165 1 : std::string refval=refstr;
166 8 : for(unsigned i=1; i<(av->copyOutput(0))->getShape()[0]; ++i) {
167 14 : refstr += "," + refval;
168 : }
169 1 : } else if( nref!=(av->copyOutput(0))->getShape()[0] ) {
170 0 : error("mismatch between number of reference values and number of ATOMS specified");
171 : }
172 2 : if( ncoeff==0 ) {
173 2 : std::string coeff=coeffstr;
174 16 : for(unsigned i=1; i<(av->copyOutput(0))->getShape()[0]; ++i) {
175 28 : coeffstr += "," + coeff;
176 : }
177 0 : } else if( ncoeff!=(av->copyOutput(0))->getShape()[0] ) {
178 0 : error("mismatch between number of coefficients and number of ATOMS specified");
179 : }
180 4 : readInputLine( getShortcutLabel() + "_ref: CONSTANT VALUES=" + refstr );
181 4 : readInputLine( getShortcutLabel() + "_coeff: CONSTANT VALUES=" + coeffstr );
182 : // Caculate difference from reference using combine
183 4 : readInputLine( getShortcutLabel() + "_comb: COMBINE ARG=" + getShortcutLabel() + "_torsions," + getShortcutLabel() + "_ref COEFFICIENTS=1,-1 PERIODIC=NO" );
184 : // Now matheval for cosine bit
185 4 : readInputLine( getShortcutLabel() + "_cos: CUSTOM ARG=" + getShortcutLabel() + "_comb," + getShortcutLabel() + "_coeff FUNC=y*(0.5+0.5*cos(x)) PERIODIC=NO");
186 : // And combine to get final value
187 4 : readInputLine( getShortcutLabel() + ": SUM ARG=" + getShortcutLabel() + "_cos PERIODIC=NO");
188 2 : }
189 :
190 : }
191 : }
|
