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