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