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