Line data Source code
1 : /* +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
2 : Copyright (c) 2013-2019 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 torsions 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 : MOLINFO MOLTYPE=protein STRUCTURE=myprotein.pdb
80 : ALPHABETA ...
81 : ATOMS1=@phi-3 REFERENCE=3.14
82 : ATOMS2=@psi-3
83 : ATOMS3=@phi-4
84 : LABEL=ab
85 : ... ALPHABETA
86 : PRINT ARG=ab FILE=colvar STRIDE=10
87 : \endplumedfile
88 :
89 : Here, \@phi-3 tells plumed that you would like to calculate the \f$\phi\f$ angle in the third residue of the protein.
90 : Similarly \@psi-4 tells plumed that you want to calculate the \f$\psi\f$ angle of the 4th residue of the protein.
91 :
92 :
93 : */
94 : //+ENDPLUMEDOC
95 :
96 3 : class AlphaBeta : public MultiColvarBase {
97 : private:
98 : std::vector<double> target;
99 : public:
100 : static void registerKeywords( Keywords& keys );
101 : explicit AlphaBeta(const ActionOptions&);
102 : virtual double compute( const unsigned& tindex, AtomValuePack& myatoms ) const ;
103 0 : bool isPeriodic() { return false; }
104 : };
105 :
106 6453 : PLUMED_REGISTER_ACTION(AlphaBeta,"ALPHABETA")
107 :
108 2 : void AlphaBeta::registerKeywords( Keywords& keys ) {
109 2 : MultiColvarBase::registerKeywords( keys );
110 8 : 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 6 : keys.reset_style("ATOMS","atoms");
116 8 : 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 torsions");
118 6 : keys.reset_style("REFERENCE","compulsory");
119 2 : }
120 :
121 1 : AlphaBeta::AlphaBeta(const ActionOptions&ao):
122 : Action(ao),
123 1 : MultiColvarBase(ao)
124 : {
125 : // Read in the atoms
126 : std::vector<AtomNumber> all_atoms;
127 2 : readAtomsLikeKeyword( "ATOMS", 4, all_atoms );
128 1 : setupMultiColvarBase( all_atoms );
129 : // Resize target
130 1 : target.resize( getFullNumberOfTasks() );
131 : // Setup central atom indices
132 1 : std::vector<bool> catom_ind(4, false);
133 : catom_ind[1]=catom_ind[2]=true;
134 1 : setAtomsForCentralAtom( catom_ind );
135 :
136 : // Read in reference values
137 : unsigned ntarget=0;
138 2 : for(unsigned i=0; i<target.size(); ++i) {
139 3 : if( !parseNumbered( "REFERENCE", i+1, target[i] ) ) break;
140 : ntarget++;
141 : }
142 1 : if( ntarget==0 ) {
143 3 : parse("REFERENCE",target[0]);
144 23 : for(unsigned i=1; i<target.size(); ++i) target[i]=target[0];
145 0 : } else if( ntarget!=target.size() ) {
146 0 : error("found wrong number of REFERENCE values");
147 : }
148 :
149 : // And setup the ActionWithVessel
150 1 : if( getNumberOfVessels()==0 ) {
151 : std::string fake_input;
152 2 : addVessel( "SUM", fake_input, -1 ); // -1 here means that this value will be named getLabel()
153 1 : readVesselKeywords(); // This makes sure resizing is done
154 : }
155 :
156 : // And check everything has been read in correctly
157 1 : checkRead();
158 1 : }
159 :
160 40 : double AlphaBeta::compute( const unsigned& tindex, AtomValuePack& myatoms ) const {
161 80 : const Vector d0=getSeparation(myatoms.getPosition(1),myatoms.getPosition(0));
162 40 : const Vector d1=getSeparation(myatoms.getPosition(2),myatoms.getPosition(1));
163 40 : const Vector d2=getSeparation(myatoms.getPosition(3),myatoms.getPosition(2));
164 :
165 40 : Vector dd0,dd1,dd2;
166 : PLMD::Torsion t;
167 40 : const double value = t.compute(d0,d1,d2,dd0,dd1,dd2);
168 80 : const double svalue = -0.5*sin(value-target[tindex]);
169 40 : const double cvalue = 1.+cos(value-target[tindex]);
170 :
171 40 : dd0 *= svalue;
172 40 : dd1 *= svalue;
173 40 : dd2 *= svalue;
174 :
175 40 : addAtomDerivatives(1, 0, dd0, myatoms);
176 40 : addAtomDerivatives(1, 1, dd1-dd0, myatoms);
177 40 : addAtomDerivatives(1, 2, dd2-dd1, myatoms);
178 40 : addAtomDerivatives(1, 3, -dd2, myatoms);
179 :
180 40 : myatoms.addBoxDerivatives(1, -(extProduct(d0,dd0)+extProduct(d1,dd1)+extProduct(d2,dd2)));
181 :
182 40 : return 0.5*cvalue;
183 : }
184 :
185 : }
186 4839 : }
|