Line data Source code
1 : /* +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
2 : Copyright (c) 2013-2020 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 "Colvar.h"
23 : #include "ColvarShortcut.h"
24 : #include "MultiColvarTemplate.h"
25 : #include "tools/Torsion.h"
26 : #include "core/ActionRegister.h"
27 :
28 : #include <string>
29 : #include <cmath>
30 :
31 : namespace PLMD {
32 : namespace colvar {
33 :
34 : //+PLUMEDOC COLVAR DIHEDRAL_CORRELATION
35 : /*
36 : Measure the correlation between a pair of dihedral angles
37 :
38 : This CV measures the correlation between two dihedral angles, $\phi$ and $\psi$, as follows:
39 :
40 : $$
41 : s = \frac{1}{2} \left[ 1 + \cos( \phi - \psi ) \right]
42 : $$
43 :
44 : An example input that computes and calculates this quantity with $\phi$ being the dihedral
45 : angle involving atoms 1, 2, 3 and 4 and $\psi$ being the angle involving atoms 7, 8, 9 and 10
46 : is shown below:
47 :
48 : ```plumed
49 : d: DIHEDRAL_CORRELATION ATOMS=1,2,3,4,5,6,7,8
50 : PRINT ARG=d FILE=colvar
51 : ```
52 :
53 : If you want to calculate the dihedral correlations between multiple pairs of dihedral angles using
54 : this action you would use an input like this one shown below:
55 :
56 : ```plumed
57 : d: DIHEDRAL_CORRELATION ...
58 : ATOMS1=1,2,3,4,5,6,7,8
59 : ATOMS2=9,10,11,12,13,14,15,16
60 : ...
61 : PRINT ARG=d FILE=colvar
62 : ```
63 :
64 : This input calculates and outputs a two dimensional vector that contains two of these dihedral correlation
65 : values. Commands similar to these are used within the [DIHCOR](DIHCOR.md) shortcut.
66 :
67 : */
68 : //+ENDPLUMEDOC
69 :
70 : //+PLUMEDOC COLVAR DIHEDRAL_CORRELATION_SCALAR
71 : /*
72 : Measure the correlation between a multiple pairs of dihedral angles
73 :
74 :
75 : \par Examples
76 :
77 : */
78 : //+ENDPLUMEDOC
79 :
80 : //+PLUMEDOC COLVAR DIHEDRAL_CORRELATION_VECTOR
81 : /*
82 : Measure the correlation between a multiple pairs of dihedral angles
83 :
84 :
85 : \par Examples
86 :
87 : */
88 : //+ENDPLUMEDOC
89 :
90 : class DihedralCorrelation : public Colvar {
91 : private:
92 : bool pbc;
93 : std::vector<double> value, masses, charges;
94 : std::vector<std::vector<Vector> > derivs;
95 : std::vector<Tensor> virial;
96 : public:
97 : static void registerKeywords( Keywords& keys );
98 : explicit DihedralCorrelation(const ActionOptions&);
99 : static void parseAtomList( const int& num, std::vector<AtomNumber>& t, ActionAtomistic* aa );
100 : static unsigned getModeAndSetupValues( ActionWithValue* av );
101 : void calculate() override;
102 : static void calculateCV( const unsigned& mode, const std::vector<double>& masses, const std::vector<double>& charges,
103 : const std::vector<Vector>& pos, std::vector<double>& vals, std::vector<std::vector<Vector> >& derivs,
104 : std::vector<Tensor>& virial, const ActionAtomistic* aa );
105 : };
106 :
107 : typedef ColvarShortcut<DihedralCorrelation> DihedralCorrelationShortcut;
108 : PLUMED_REGISTER_ACTION(DihedralCorrelationShortcut,"DIHEDRAL_CORRELATION")
109 : PLUMED_REGISTER_ACTION(DihedralCorrelation,"DIHEDRAL_CORRELATION_SCALAR")
110 : typedef MultiColvarTemplate<DihedralCorrelation> DihedralCorrelationMulti;
111 : PLUMED_REGISTER_ACTION(DihedralCorrelationMulti,"DIHEDRAL_CORRELATION_VECTOR")
112 :
113 10 : void DihedralCorrelation::registerKeywords( Keywords& keys ) {
114 10 : Colvar::registerKeywords( keys );
115 20 : keys.setDisplayName("DIHEDRAL_CORRELATION");
116 10 : keys.add("atoms","ATOMS","the set of 8 atoms that are being used to calculate this quantity");
117 10 : keys.add("hidden","NO_ACTION_LOG","suppresses printing from action on the log");
118 20 : keys.setValueDescription("scalar/vector","the DIHEDRAL_CORRELATION for these atoms");
119 10 : }
120 :
121 0 : DihedralCorrelation::DihedralCorrelation(const ActionOptions&ao):
122 : PLUMED_COLVAR_INIT(ao),
123 0 : pbc(true),
124 0 : value(1),
125 0 : derivs(1),
126 0 : virial(1) {
127 0 : derivs[0].resize(8);
128 : std::vector<AtomNumber> atoms;
129 0 : parseAtomList(-1,atoms,this);
130 0 : if( atoms.size()!=8 ) {
131 0 : error("Number of specified atoms should be 8");
132 : }
133 :
134 0 : bool nopbc=!pbc;
135 0 : parseFlag("NOPBC",nopbc);
136 0 : pbc=!nopbc;
137 :
138 0 : if(pbc) {
139 0 : log.printf(" using periodic boundary conditions\n");
140 : } else {
141 0 : log.printf(" without periodic boundary conditions\n");
142 : }
143 0 : addValueWithDerivatives();
144 0 : setNotPeriodic();
145 0 : }
146 :
147 3 : void DihedralCorrelation::parseAtomList( const int& num, std::vector<AtomNumber>& t, ActionAtomistic* aa ) {
148 3 : aa->parseAtomList("ATOMS",num,t);
149 3 : if( num<0 && t.size()!=8 ) {
150 0 : aa->error("Number of specified atoms should be 8");
151 : }
152 3 : if( t.size()==8 ) {
153 2 : aa->log.printf(" correlation between dihedral angle for atoms %d %d %d %d and atoms %d %d %d %d\n",
154 : t[0].serial(),t[1].serial(),t[2].serial(),t[3].serial(),t[4].serial(),t[5].serial(),t[6].serial(),t[7].serial());
155 : }
156 3 : }
157 :
158 1 : unsigned DihedralCorrelation::getModeAndSetupValues( ActionWithValue* av ) {
159 1 : av->addValueWithDerivatives();
160 1 : av->setNotPeriodic();
161 1 : return 0;
162 : }
163 :
164 0 : void DihedralCorrelation::calculate() {
165 :
166 0 : if(pbc) {
167 0 : makeWhole();
168 : }
169 0 : calculateCV( 0, masses, charges, getPositions(), value, derivs, virial, this );
170 0 : setValue( value[0] );
171 0 : for(unsigned i=0; i<derivs[0].size(); ++i) {
172 0 : setAtomsDerivatives( i, derivs[0][i] );
173 : }
174 0 : setBoxDerivatives( virial[0] );
175 0 : }
176 :
177 10 : void DihedralCorrelation::calculateCV( const unsigned& mode, const std::vector<double>& masses, const std::vector<double>& charges,
178 : const std::vector<Vector>& pos, std::vector<double>& vals, std::vector<std::vector<Vector> >& derivs,
179 : std::vector<Tensor>& virial, const ActionAtomistic* aa ) {
180 10 : const Vector d10=delta(pos[1],pos[0]);
181 10 : const Vector d11=delta(pos[2],pos[1]);
182 10 : const Vector d12=delta(pos[3],pos[2]);
183 :
184 10 : Vector dd10,dd11,dd12;
185 : PLMD::Torsion t1;
186 10 : const double phi1 = t1.compute(d10,d11,d12,dd10,dd11,dd12);
187 :
188 10 : const Vector d20=delta(pos[5],pos[4]);
189 10 : const Vector d21=delta(pos[6],pos[5]);
190 10 : const Vector d22=delta(pos[7],pos[6]);
191 :
192 10 : Vector dd20,dd21,dd22;
193 : PLMD::Torsion t2;
194 10 : const double phi2 = t2.compute( d20, d21, d22, dd20, dd21, dd22 );
195 :
196 : // Calculate value
197 10 : const double diff = phi2 - phi1;
198 10 : vals[0] = 0.5*(1.+std::cos(diff));
199 : // Derivatives wrt phi1
200 10 : const double dval = 0.5*std::sin(diff);
201 10 : dd10 *= dval;
202 10 : dd11 *= dval;
203 10 : dd12 *= dval;
204 : // And add
205 10 : derivs[0][0]=dd10;
206 10 : derivs[0][1]=dd11-dd10;
207 10 : derivs[0][2]=dd12-dd11;
208 10 : derivs[0][3]=-dd12;
209 : // Derivative wrt phi2
210 10 : dd20 *= -dval;
211 10 : dd21 *= -dval;
212 10 : dd22 *= -dval;
213 : // And add
214 10 : derivs[0][4]=dd20;
215 10 : derivs[0][5]=dd21-dd20;
216 10 : derivs[0][6]=dd22-dd21;
217 10 : derivs[0][7]=-dd22;
218 10 : virial[0] = -(extProduct(d10,dd10)+extProduct(d11,dd11)+extProduct(d12,dd12)) - (extProduct(d20,dd20)+extProduct(d21,dd21)+extProduct(d22,dd22));
219 10 : }
220 :
221 : }
222 : }
|
