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 {
104 0 : return false;
105 : }
106 : };
107 :
108 13787 : PLUMED_REGISTER_ACTION(AlphaBeta,"ALPHABETA")
109 :
110 5 : void AlphaBeta::registerKeywords( Keywords& keys ) {
111 5 : MultiColvarBase::registerKeywords( keys );
112 10 : keys.add("numbered","ATOMS","the atoms involved in each of the alpha-beta variables you wish to calculate. "
113 : "Keywords like ATOMS1, ATOMS2, ATOMS3,... should be listed and one alpha-beta values will be "
114 : "calculated for each ATOM keyword you specify (all ATOM keywords should "
115 : "specify the indices of four atoms). The eventual number of quantities calculated by this "
116 : "action will depend on what functions of the distribution you choose to calculate.");
117 10 : keys.reset_style("ATOMS","atoms");
118 10 : keys.add("numbered","REFERENCE","the reference values for each of the torsional angles. If you use a single REFERENCE value the "
119 : "same reference value is used for all torsional angles");
120 10 : keys.add("numbered","COEFFICIENT","the coefficient for each of the torsional angles. If you use a single COEFFICIENT value the "
121 : "same reference value is used for all torsional angles");
122 10 : keys.reset_style("REFERENCE","compulsory");
123 10 : keys.reset_style("COEFFICIENT","optional");
124 5 : }
125 :
126 1 : AlphaBeta::AlphaBeta(const ActionOptions&ao):
127 : Action(ao),
128 1 : MultiColvarBase(ao) {
129 : // Read in the atoms
130 : std::vector<AtomNumber> all_atoms;
131 1 : readAtomsLikeKeyword( "ATOMS", 4, all_atoms );
132 1 : setupMultiColvarBase( all_atoms );
133 : // Resize target
134 1 : target.resize( getFullNumberOfTasks() );
135 : // Resize coeff
136 1 : coefficient.resize( getFullNumberOfTasks(), 1.0);
137 : // Setup central atom indices
138 1 : std::vector<bool> catom_ind(4, false);
139 1 : catom_ind[1]=catom_ind[2]=true;
140 1 : setAtomsForCentralAtom( catom_ind );
141 :
142 : // Read in reference values
143 : unsigned ntarget=0;
144 1 : for(unsigned i=0; i<target.size(); ++i) {
145 2 : if( !parseNumbered( "REFERENCE", i+1, target[i] ) ) {
146 : break;
147 : }
148 0 : ntarget++;
149 : }
150 1 : if( ntarget==0 ) {
151 1 : parse("REFERENCE",target[0]);
152 8 : for(unsigned i=1; i<target.size(); ++i) {
153 7 : target[i]=target[0];
154 : }
155 0 : } else if( ntarget!=target.size() ) {
156 0 : error("found wrong number of REFERENCE values");
157 : }
158 :
159 : // Read in reference values
160 : unsigned ncoefficient=0;
161 1 : for(unsigned i=0; i<coefficient.size(); ++i) {
162 2 : if( !parseNumbered( "COEFFICIENT", i+1, coefficient[i] ) ) {
163 : break;
164 : }
165 0 : ncoefficient++;
166 : }
167 1 : if( ncoefficient==0 ) {
168 1 : parse("COEFFICIENT",coefficient[0]);
169 8 : for(unsigned i=1; i<coefficient.size(); ++i) {
170 7 : coefficient[i]=coefficient[0];
171 : }
172 0 : } else if( ncoefficient !=coefficient.size() ) {
173 0 : error("found wrong number of COEFFICIENT values");
174 : }
175 :
176 : // And setup the ActionWithVessel
177 1 : if( getNumberOfVessels()==0 ) {
178 : std::string fake_input;
179 1 : addVessel( "SUM", fake_input, -1 ); // -1 here means that this value will be named getLabel()
180 1 : readVesselKeywords(); // This makes sure resizing is done
181 : }
182 :
183 : // And check everything has been read in correctly
184 1 : checkRead();
185 1 : }
186 :
187 40 : double AlphaBeta::compute( const unsigned& tindex, AtomValuePack& myatoms ) const {
188 40 : const Vector d0=getSeparation(myatoms.getPosition(1),myatoms.getPosition(0));
189 40 : const Vector d1=getSeparation(myatoms.getPosition(2),myatoms.getPosition(1));
190 40 : const Vector d2=getSeparation(myatoms.getPosition(3),myatoms.getPosition(2));
191 :
192 40 : Vector dd0,dd1,dd2;
193 : PLMD::Torsion t;
194 40 : const double value = t.compute(d0,d1,d2,dd0,dd1,dd2);
195 40 : const double svalue = -0.5*coefficient[tindex]*std::sin(value-target[tindex]);
196 40 : const double cvalue = coefficient[tindex]*(1.+std::cos(value-target[tindex]));
197 :
198 40 : dd0 *= svalue;
199 40 : dd1 *= svalue;
200 40 : dd2 *= svalue;
201 :
202 40 : addAtomDerivatives(1, 0, dd0, myatoms);
203 40 : addAtomDerivatives(1, 1, dd1-dd0, myatoms);
204 40 : addAtomDerivatives(1, 2, dd2-dd1, myatoms);
205 40 : addAtomDerivatives(1, 3, -dd2, myatoms);
206 :
207 40 : myatoms.addBoxDerivatives(1, -(extProduct(d0,dd0)+extProduct(d1,dd1)+extProduct(d2,dd2)));
208 :
209 40 : return 0.5*cvalue;
210 : }
211 :
212 : }
213 : }
|
