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 DIHCOR
36 : /*
37 : Measures the degree of similarity between dihedral angles.
38 :
39 : This colvar calculates the following quantity.
40 :
41 : \f[
42 : s = \frac{1}{2} \sum_i \left[ 1 + \cos( \phi_i - \psi_i ) \right]
43 : \f]
44 :
45 : where the \f$\phi_i\f$ and \f$\psi\f$ values and the instantaneous values for the \ref TORSION angles of interest.
46 :
47 : \par Examples
48 :
49 : The following provides an example input for the dihcor action
50 :
51 : \plumedfile
52 : DIHCOR ...
53 : ATOMS1=1,2,3,4,5,6,7,8
54 : ATOMS2=5,6,7,8,9,10,11,12
55 : LABEL=dih
56 : ... DIHCOR
57 : PRINT ARG=dih FILE=colvar STRIDE=10
58 : \endplumedfile
59 :
60 : In the above input we are calculating the correation between the torsion angle involving atoms 1, 2, 3 and 4 and the torsion angle
61 : involving atoms 5, 6, 7 and 8. This is then added to the correlation betwene the torsion angle involving atoms 5, 6, 7 and 8 and the
62 : correlation angle involving atoms 9, 10, 11 and 12.
63 :
64 : Writing out the atoms involved in all the torsions in this way can be rather tedious. Thankfully if you are working with protein you
65 : can avoid this by using the \ref MOLINFO command. PLUMED uses the pdb file that you provide to this command to learn
66 : about the topology of the protein molecule. This means that you can specify torsion angles using the following syntax:
67 :
68 : \plumedfile
69 : MOLINFO MOLTYPE=protein STRUCTURE=myprotein.pdb
70 : DIHCOR ...
71 : ATOMS1=@phi-3,@psi-3
72 : ATOMS2=@psi-3,@phi-4
73 : ATOMS4=@phi-4,@psi-4
74 : ... DIHCOR
75 : PRINT ARG=dih FILE=colvar STRIDE=10
76 : \endplumedfile
77 :
78 : Here, \@phi-3 tells plumed that you would like to calculate the \f$\phi\f$ angle in the third residue of the protein.
79 : Similarly \@psi-4 tells plumed that you want to calculate the \f$\psi\f$ angle of the 4th residue of the protein.
80 :
81 : */
82 : //+ENDPLUMEDOC
83 :
84 4 : class DihedralCorrelation : public MultiColvarBase {
85 : private:
86 : public:
87 : static void registerKeywords( Keywords& keys );
88 : explicit DihedralCorrelation(const ActionOptions&);
89 : virtual double compute( const unsigned& tindex, AtomValuePack& myatoms ) const ;
90 0 : bool isPeriodic() { return false; }
91 : };
92 :
93 6454 : PLUMED_REGISTER_ACTION(DihedralCorrelation,"DIHCOR")
94 :
95 3 : void DihedralCorrelation::registerKeywords( Keywords& keys ) {
96 3 : MultiColvarBase::registerKeywords( keys );
97 12 : keys.add("numbered","ATOMS","the atoms involved in each of the dihedral correlation values you wish to calculate. "
98 : "Keywords like ATOMS1, ATOMS2, ATOMS3,... should be listed and one dihedral correlation will be "
99 : "calculated for each ATOM keyword you specify (all ATOM keywords should "
100 : "specify the indices of 8 atoms). The eventual number of quantities calculated by this "
101 : "action will depend on what functions of the distribution you choose to calculate.");
102 9 : keys.reset_style("ATOMS","atoms");
103 3 : }
104 :
105 2 : DihedralCorrelation::DihedralCorrelation(const ActionOptions&ao):
106 : Action(ao),
107 2 : MultiColvarBase(ao)
108 : {
109 : // Read in the atoms
110 : std::vector<AtomNumber> all_atoms;
111 4 : readAtomsLikeKeyword( "ATOMS", 8, all_atoms );
112 2 : setupMultiColvarBase( all_atoms );
113 : // Stuff for central atoms
114 2 : std::vector<bool> catom_ind(8, false);
115 : catom_ind[1]=catom_ind[2]=catom_ind[5]=catom_ind[6]=true;
116 2 : setAtomsForCentralAtom( catom_ind );
117 :
118 : // And setup the ActionWithVessel
119 2 : if( getNumberOfVessels()==0 ) {
120 : std::string fake_input;
121 4 : addVessel( "SUM", fake_input, -1 ); // -1 here means that this value will be named getLabel()
122 2 : readVesselKeywords(); // This makes sure resizing is done
123 : }
124 :
125 : // And check everything has been read in correctly
126 2 : checkRead();
127 2 : }
128 :
129 590 : double DihedralCorrelation::compute( const unsigned& tindex, AtomValuePack& myatoms ) const {
130 1180 : const Vector d10=getSeparation(myatoms.getPosition(1),myatoms.getPosition(0));
131 590 : const Vector d11=getSeparation(myatoms.getPosition(2),myatoms.getPosition(1));
132 590 : const Vector d12=getSeparation(myatoms.getPosition(3),myatoms.getPosition(2));
133 :
134 590 : Vector dd10,dd11,dd12;
135 : PLMD::Torsion t1;
136 590 : const double phi1 = t1.compute(d10,d11,d12,dd10,dd11,dd12);
137 :
138 590 : const Vector d20=getSeparation(myatoms.getPosition(5),myatoms.getPosition(4));
139 590 : const Vector d21=getSeparation(myatoms.getPosition(6),myatoms.getPosition(5));
140 590 : const Vector d22=getSeparation(myatoms.getPosition(7),myatoms.getPosition(6));
141 :
142 590 : Vector dd20,dd21,dd22;
143 : PLMD::Torsion t2;
144 590 : const double phi2 = t2.compute( d20, d21, d22, dd20, dd21, dd22 );
145 :
146 : // Calculate value
147 590 : const double diff = phi2 - phi1;
148 590 : const double value = 0.5*(1.+cos(diff));
149 : // Derivatives wrt phi1
150 590 : const double dval = 0.5*sin(diff);
151 590 : dd10 *= dval;
152 590 : dd11 *= dval;
153 590 : dd12 *= dval;
154 : // And add
155 590 : addAtomDerivatives(1, 0, dd10, myatoms );
156 590 : addAtomDerivatives(1, 1, dd11-dd10, myatoms );
157 590 : addAtomDerivatives(1, 2, dd12-dd11, myatoms );
158 590 : addAtomDerivatives(1, 3, -dd12, myatoms );
159 590 : myatoms.addBoxDerivatives (1, -(extProduct(d10,dd10)+extProduct(d11,dd11)+extProduct(d12,dd12)));
160 : // Derivative wrt phi2
161 590 : dd20 *= -dval;
162 590 : dd21 *= -dval;
163 590 : dd22 *= -dval;
164 : // And add
165 590 : addAtomDerivatives(1, 4, dd20, myatoms );
166 590 : addAtomDerivatives(1, 5, dd21-dd20, myatoms );
167 590 : addAtomDerivatives(1, 6, dd22-dd21, myatoms );
168 590 : addAtomDerivatives(1, 7, -dd22, myatoms );
169 590 : myatoms.addBoxDerivatives(1, -(extProduct(d20,dd20)+extProduct(d21,dd21)+extProduct(d22,dd22)));
170 :
171 590 : return value;
172 : }
173 :
174 : }
175 4839 : }
|
