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 "MultiColvarBase.h"
23 : #include "AtomValuePack.h"
24 : #include "tools/Torsion.h"
25 : #include "core/ActionRegister.h"
26 :
27 : #include <string>
28 : #include <cmath>
29 :
30 : namespace PLMD {
31 : namespace multicolvar {
32 :
33 : //+PLUMEDOC COLVAR ALPHABETA
34 : /*
35 : Measures a distance including pbc between the instantaneous values of a set of torsional angles and set of reference values.
36 :
37 : This colvar calculates the following quantity.
38 :
39 : \f[
40 : s = \frac{1}{2} \sum_i \left[ 1 + \cos( \phi_i - \phi_i^{\textrm{Ref}} ) \right]
41 : \f]
42 :
43 : where the \f$\phi_i\f$ values are the instantaneous values for the \ref TORSION angles of interest.
44 : The \f$\phi_i^{\textrm{Ref}}\f$ values are the user-specified reference values for the torsional angles.
45 :
46 : \par Examples
47 :
48 : The following provides an example of the input for an alpha beta similarity.
49 :
50 : \plumedfile
51 : ALPHABETA ...
52 : ATOMS1=168,170,172,188 REFERENCE1=3.14
53 : ATOMS2=170,172,188,190 REFERENCE2=3.14
54 : ATOMS3=188,190,192,230 REFERENCE3=3.14
55 : LABEL=ab
56 : ... ALPHABETA
57 : PRINT ARG=ab FILE=colvar STRIDE=10
58 : \endplumedfile
59 :
60 : Because all the reference values are the same we can calculate the same quantity using
61 :
62 : \plumedfile
63 : ALPHABETA ...
64 : ATOMS1=168,170,172,188 REFERENCE=3.14
65 : ATOMS2=170,172,188,190
66 : ATOMS3=188,190,192,230
67 : LABEL=ab
68 : ... ALPHABETA
69 : PRINT ARG=ab FILE=colvar STRIDE=10
70 : \endplumedfile
71 :
72 : 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
73 : can avoid this by using the \ref MOLINFO command. PLUMED uses the pdb file that you provide to this command to learn
74 : about the topology of the protein molecule. This means that you can specify torsion angles using the following syntax:
75 :
76 : \plumedfile
77 : #SETTINGS MOLFILE=regtest/basic/rt32/helix.pdb
78 : MOLINFO MOLTYPE=protein STRUCTURE=myprotein.pdb
79 : ALPHABETA ...
80 : ATOMS1=@phi-3 REFERENCE=3.14
81 : ATOMS2=@psi-3
82 : ATOMS3=@phi-4
83 : LABEL=ab
84 : ... ALPHABETA
85 : PRINT ARG=ab FILE=colvar STRIDE=10
86 : \endplumedfile
87 :
88 : Here, \@phi-3 tells plumed that you would like to calculate the \f$\phi\f$ angle in the third residue of the protein.
89 : Similarly \@psi-4 tells plumed that you want to calculate the \f$\psi\f$ angle of the fourth residue of the protein.
90 :
91 :
92 : */
93 : //+ENDPLUMEDOC
94 :
95 : class AlphaBeta : public MultiColvarBase {
96 : private:
97 : std::vector<double> target;
98 : std::vector<double> coefficient;
99 : public:
100 : static void registerKeywords( Keywords& keys );
101 : explicit AlphaBeta(const ActionOptions&);
102 : double compute( const unsigned& tindex, AtomValuePack& myatoms ) const override;
103 0 : bool isPeriodic() override { return false; }
104 : };
105 :
106 10421 : PLUMED_REGISTER_ACTION(AlphaBeta,"ALPHABETA")
107 :
108 2 : void AlphaBeta::registerKeywords( Keywords& keys ) {
109 2 : MultiColvarBase::registerKeywords( keys );
110 4 : keys.add("numbered","ATOMS","the atoms involved in each of the alpha-beta variables you wish to calculate. "
111 : "Keywords like ATOMS1, ATOMS2, ATOMS3,... should be listed and one alpha-beta values will be "
112 : "calculated for each ATOM keyword you specify (all ATOM keywords should "
113 : "specify the indices of four atoms). The eventual number of quantities calculated by this "
114 : "action will depend on what functions of the distribution you choose to calculate.");
115 4 : keys.reset_style("ATOMS","atoms");
116 4 : keys.add("numbered","REFERENCE","the reference values for each of the torsional angles. If you use a single REFERENCE value the "
117 : "same reference value is used for all torsional angles");
118 4 : keys.add("numbered","COEFFICIENT","the coefficient for each of the torsional angles. If you use a single COEFFICIENT value the "
119 : "same reference value is used for all torsional angles");
120 4 : keys.reset_style("REFERENCE","compulsory");
121 4 : keys.reset_style("COEFFICIENT","optional");
122 2 : }
123 :
124 1 : AlphaBeta::AlphaBeta(const ActionOptions&ao):
125 : Action(ao),
126 1 : MultiColvarBase(ao)
127 : {
128 : // Read in the atoms
129 : std::vector<AtomNumber> all_atoms;
130 1 : readAtomsLikeKeyword( "ATOMS", 4, all_atoms );
131 1 : setupMultiColvarBase( all_atoms );
132 : // Resize target
133 1 : target.resize( getFullNumberOfTasks() );
134 : // Resize coeff
135 1 : coefficient.resize( getFullNumberOfTasks(), 1.0);
136 : // Setup central atom indices
137 1 : std::vector<bool> catom_ind(4, false);
138 1 : catom_ind[1]=catom_ind[2]=true;
139 1 : setAtomsForCentralAtom( catom_ind );
140 :
141 : // Read in reference values
142 : unsigned ntarget=0;
143 1 : for(unsigned i=0; i<target.size(); ++i) {
144 2 : if( !parseNumbered( "REFERENCE", i+1, target[i] ) ) break;
145 0 : ntarget++;
146 : }
147 1 : if( ntarget==0 ) {
148 1 : parse("REFERENCE",target[0]);
149 8 : for(unsigned i=1; i<target.size(); ++i) target[i]=target[0];
150 0 : } else if( ntarget!=target.size() ) {
151 0 : error("found wrong number of REFERENCE values");
152 : }
153 :
154 : // Read in reference values
155 : unsigned ncoefficient=0;
156 1 : for(unsigned i=0; i<coefficient.size(); ++i) {
157 2 : if( !parseNumbered( "COEFFICIENT", i+1, coefficient[i] ) ) break;
158 0 : ncoefficient++;
159 : }
160 1 : if( ncoefficient==0 ) {
161 1 : parse("COEFFICIENT",coefficient[0]);
162 8 : for(unsigned i=1; i<coefficient.size(); ++i) coefficient[i]=coefficient[0];
163 0 : } else if( ncoefficient !=coefficient.size() ) {
164 0 : error("found wrong number of COEFFICIENT values");
165 : }
166 :
167 : // And setup the ActionWithVessel
168 1 : if( getNumberOfVessels()==0 ) {
169 : std::string fake_input;
170 1 : addVessel( "SUM", fake_input, -1 ); // -1 here means that this value will be named getLabel()
171 1 : readVesselKeywords(); // This makes sure resizing is done
172 : }
173 :
174 : // And check everything has been read in correctly
175 1 : checkRead();
176 1 : }
177 :
178 40 : double AlphaBeta::compute( const unsigned& tindex, AtomValuePack& myatoms ) const {
179 40 : const Vector d0=getSeparation(myatoms.getPosition(1),myatoms.getPosition(0));
180 40 : const Vector d1=getSeparation(myatoms.getPosition(2),myatoms.getPosition(1));
181 40 : const Vector d2=getSeparation(myatoms.getPosition(3),myatoms.getPosition(2));
182 :
183 40 : Vector dd0,dd1,dd2;
184 : PLMD::Torsion t;
185 40 : const double value = t.compute(d0,d1,d2,dd0,dd1,dd2);
186 40 : const double svalue = -0.5*coefficient[tindex]*std::sin(value-target[tindex]);
187 40 : const double cvalue = coefficient[tindex]*(1.+std::cos(value-target[tindex]));
188 :
189 40 : dd0 *= svalue;
190 40 : dd1 *= svalue;
191 40 : dd2 *= svalue;
192 :
193 40 : addAtomDerivatives(1, 0, dd0, myatoms);
194 40 : addAtomDerivatives(1, 1, dd1-dd0, myatoms);
195 40 : addAtomDerivatives(1, 2, dd2-dd1, myatoms);
196 40 : addAtomDerivatives(1, 3, -dd2, myatoms);
197 :
198 40 : myatoms.addBoxDerivatives(1, -(extProduct(d0,dd0)+extProduct(d1,dd1)+extProduct(d2,dd2)));
199 :
200 40 : return 0.5*cvalue;
201 : }
202 :
203 : }
204 : }
|
