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