Line data Source code
1 : /* +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
2 : Copyright (c) 2013-2018 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 "MultiColvar.h"
23 : #include "tools/Torsion.h"
24 : #include "core/ActionRegister.h"
25 :
26 : #include <string>
27 : #include <cmath>
28 :
29 : using namespace std;
30 :
31 : namespace PLMD {
32 : namespace multicolvar {
33 :
34 : //+PLUMEDOC COLVAR ALPHABETA
35 : /*
36 : Measures a distance including pbc between the instantaneous values of a set of torsional angles and set of reference values.
37 :
38 : This colvar calculates the following quantity.
39 :
40 : \f[
41 : s = \frac{1}{2} \sum_i \left[ 1 + \cos( \phi_i - \phi_i^{\textrm{Ref}} ) \right]
42 : \f]
43 :
44 : where the \f$\phi_i\f$ values are the instantaneous values for the \ref TORSION angles of interest.
45 : The \f$\phi_i^{\textrm{Ref}}\f$ values are the user-specified reference values for the torsional angles.
46 :
47 : \par Examples
48 :
49 : The following provides an example of the input for an alpha beta similarity.
50 :
51 : \verbatim
52 : ALPHABETA ...
53 : ATOMS1=168,170,172,188 REFERENCE1=3.14
54 : ATOMS2=170,172,188,190 REFERENCE2=3.14
55 : ATOMS3=188,190,192,230 REFERENCE3=3.14
56 : LABEL=ab
57 : ... ALPHABETA
58 : PRINT ARG=ab FILE=colvar STRIDE=10
59 : \endverbatim
60 :
61 : Because all the reference values are the same we can calculate the same quantity using
62 :
63 : \verbatim
64 : ALPHABETA ...
65 : ATOMS1=168,170,172,188 REFERENCE=3.14
66 : ATOMS2=170,172,188,190
67 : ATOMS3=188,190,192,230
68 : LABEL=ab
69 : ... ALPHABETA
70 : PRINT ARG=ab FILE=colvar STRIDE=10
71 : \endverbatim
72 :
73 : Writing out the atoms involved in all the torsions in this way can be rather tedious. Thankfully if you are working with protein you
74 : can avoid this by using the \ref MOLINFO command. PLUMED uses the pdb file that you provide to this command to learn
75 : about the topology of the protein molecule. This means that you can specify torsion angles using the following syntax:
76 :
77 : \verbatim
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 : \endverbatim
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 4th residue of the protein.
90 :
91 :
92 : */
93 : //+ENDPLUMEDOC
94 :
95 4 : class AlphaBeta : public MultiColvar {
96 : private:
97 : std::vector<double> target;
98 : public:
99 : static void registerKeywords( Keywords& keys );
100 : explicit AlphaBeta(const ActionOptions&);
101 : virtual double compute( const unsigned& tindex, AtomValuePack& myatoms ) const ;
102 0 : bool isPeriodic() { return false; }
103 : };
104 :
105 2525 : PLUMED_REGISTER_ACTION(AlphaBeta,"ALPHABETA")
106 :
107 3 : void AlphaBeta::registerKeywords( Keywords& keys ) {
108 3 : MultiColvar::registerKeywords( keys );
109 3 : keys.use("ATOMS");
110 : keys.add("numbered","REFERENCE","the reference values for each of the torsional angles. If you use a single REFERENCE value the "
111 3 : "same reference value is used for all torsions");
112 3 : keys.reset_style("REFERENCE","compulsory");
113 3 : }
114 :
115 2 : AlphaBeta::AlphaBeta(const ActionOptions&ao):
116 2 : PLUMED_MULTICOLVAR_INIT(ao)
117 : {
118 : // Read in the atoms
119 2 : int natoms=4; std::vector<AtomNumber> all_atoms;
120 2 : readAtoms( natoms, all_atoms );
121 : // Resize target
122 2 : target.resize( getFullNumberOfTasks() );
123 : // Setup central atom indices
124 4 : std::vector<bool> catom_ind(4, false);
125 2 : catom_ind[1]=catom_ind[2]=true;
126 2 : setAtomsForCentralAtom( catom_ind );
127 :
128 : // Read in reference values
129 2 : unsigned ntarget=0;
130 2 : for(unsigned i=0; i<target.size(); ++i) {
131 2 : if( !parseNumbered( "REFERENCE", i+1, target[i] ) ) break;
132 0 : ntarget++;
133 : }
134 2 : if( ntarget==0 ) {
135 2 : parse("REFERENCE",target[0]);
136 2 : for(unsigned i=1; i<target.size(); ++i) target[i]=target[0];
137 0 : } else if( ntarget!=target.size() ) {
138 0 : error("found wrong number of REFERENCE values");
139 : }
140 :
141 : // And setup the ActionWithVessel
142 2 : if( getNumberOfVessels()==0 ) {
143 2 : std::string fake_input;
144 2 : addVessel( "SUM", fake_input, -1 ); // -1 here means that this value will be named getLabel()
145 2 : readVesselKeywords(); // This makes sure resizing is done
146 : }
147 :
148 : // And check everything has been read in correctly
149 4 : checkRead();
150 2 : }
151 :
152 350 : double AlphaBeta::compute( const unsigned& tindex, AtomValuePack& myatoms ) const {
153 350 : const Vector d0=getSeparation(myatoms.getPosition(1),myatoms.getPosition(0));
154 350 : const Vector d1=getSeparation(myatoms.getPosition(2),myatoms.getPosition(1));
155 350 : const Vector d2=getSeparation(myatoms.getPosition(3),myatoms.getPosition(2));
156 :
157 350 : Vector dd0,dd1,dd2;
158 : PLMD::Torsion t;
159 350 : const double value = t.compute(d0,d1,d2,dd0,dd1,dd2);
160 350 : const double svalue = -0.5*sin(value-target[tindex]);
161 350 : const double cvalue = 1.+cos(value-target[tindex]);
162 :
163 350 : dd0 *= svalue;
164 350 : dd1 *= svalue;
165 350 : dd2 *= svalue;
166 :
167 350 : addAtomDerivatives(1, 0, dd0, myatoms);
168 350 : addAtomDerivatives(1, 1, dd1-dd0, myatoms);
169 350 : addAtomDerivatives(1, 2, dd2-dd1, myatoms);
170 350 : addAtomDerivatives(1, 3, -dd2, myatoms);
171 :
172 350 : myatoms.addBoxDerivatives(1, -(extProduct(d0,dd0)+extProduct(d1,dd1)+extProduct(d2,dd2)));
173 :
174 350 : return 0.5*cvalue;
175 : }
176 :
177 : }
178 2523 : }
|
