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 :
39 : \par Examples
40 :
41 : */
42 : //+ENDPLUMEDOC
43 :
44 : //+PLUMEDOC COLVAR DIHEDRAL_CORRELATION_SCALAR
45 : /*
46 : Measure the correlation between a multiple pairs of dihedral angles
47 :
48 :
49 : \par Examples
50 :
51 : */
52 : //+ENDPLUMEDOC
53 :
54 : //+PLUMEDOC COLVAR DIHEDRAL_CORRELATION_VECTOR
55 : /*
56 : Measure the correlation between a multiple pairs of dihedral angles
57 :
58 :
59 : \par Examples
60 :
61 : */
62 : //+ENDPLUMEDOC
63 :
64 : class DihedralCorrelation : public Colvar {
65 : private:
66 : bool pbc;
67 : std::vector<double> value, masses, charges;
68 : std::vector<std::vector<Vector> > derivs;
69 : std::vector<Tensor> virial;
70 : public:
71 : static void registerKeywords( Keywords& keys );
72 : explicit DihedralCorrelation(const ActionOptions&);
73 : static void parseAtomList( const int& num, std::vector<AtomNumber>& t, ActionAtomistic* aa );
74 : static unsigned getModeAndSetupValues( ActionWithValue* av );
75 : void calculate() override;
76 : static void calculateCV( const unsigned& mode, const std::vector<double>& masses, const std::vector<double>& charges,
77 : const std::vector<Vector>& pos, std::vector<double>& vals, std::vector<std::vector<Vector> >& derivs,
78 : std::vector<Tensor>& virial, const ActionAtomistic* aa );
79 : };
80 :
81 : typedef ColvarShortcut<DihedralCorrelation> DihedralCorrelationShortcut;
82 : PLUMED_REGISTER_ACTION(DihedralCorrelationShortcut,"DIHEDRAL_CORRELATION")
83 : PLUMED_REGISTER_ACTION(DihedralCorrelation,"DIHEDRAL_CORRELATION_SCALAR")
84 : typedef MultiColvarTemplate<DihedralCorrelation> DihedralCorrelationMulti;
85 : PLUMED_REGISTER_ACTION(DihedralCorrelationMulti,"DIHEDRAL_CORRELATION_VECTOR")
86 :
87 10 : void DihedralCorrelation::registerKeywords( Keywords& keys ) {
88 10 : Colvar::registerKeywords( keys );
89 10 : keys.setDisplayName("DIHEDRAL_CORRELATION");
90 20 : keys.add("atoms","ATOMS","the set of 8 atoms that are being used to calculate this quantity");
91 20 : keys.add("hidden","NO_ACTION_LOG","suppresses printing from action on the log");
92 10 : keys.setValueDescription("the DIHEDRAL_CORRELATION for these atoms");
93 10 : }
94 :
95 0 : DihedralCorrelation::DihedralCorrelation(const ActionOptions&ao):
96 : PLUMED_COLVAR_INIT(ao),
97 0 : pbc(true),
98 0 : value(1),
99 0 : derivs(1),
100 0 : virial(1) {
101 0 : derivs[0].resize(8);
102 : std::vector<AtomNumber> atoms;
103 0 : parseAtomList(-1,atoms,this);
104 0 : if( atoms.size()!=8 ) {
105 0 : error("Number of specified atoms should be 8");
106 : }
107 :
108 0 : bool nopbc=!pbc;
109 0 : parseFlag("NOPBC",nopbc);
110 0 : pbc=!nopbc;
111 :
112 0 : if(pbc) {
113 0 : log.printf(" using periodic boundary conditions\n");
114 : } else {
115 0 : log.printf(" without periodic boundary conditions\n");
116 : }
117 0 : }
118 :
119 3 : void DihedralCorrelation::parseAtomList( const int& num, std::vector<AtomNumber>& t, ActionAtomistic* aa ) {
120 3 : aa->parseAtomList("ATOMS",num,t);
121 3 : if( num<0 && t.size()!=8 ) {
122 0 : aa->error("Number of specified atoms should be 8");
123 : }
124 3 : if( t.size()==8 ) {
125 2 : aa->log.printf(" correlation between dihedral angle for atoms %d %d %d %d and atoms %d %d %d %d\n",
126 : t[0].serial(),t[1].serial(),t[2].serial(),t[3].serial(),t[4].serial(),t[5].serial(),t[6].serial(),t[7].serial());
127 : }
128 3 : }
129 :
130 1 : unsigned DihedralCorrelation::getModeAndSetupValues( ActionWithValue* av ) {
131 1 : av->addValueWithDerivatives();
132 1 : av->setNotPeriodic();
133 1 : return 0;
134 : }
135 :
136 0 : void DihedralCorrelation::calculate() {
137 :
138 0 : if(pbc) {
139 0 : makeWhole();
140 : }
141 0 : calculateCV( 0, masses, charges, getPositions(), value, derivs, virial, this );
142 0 : setValue( value[0] );
143 0 : for(unsigned i=0; i<derivs[0].size(); ++i) {
144 0 : setAtomsDerivatives( i, derivs[0][i] );
145 : }
146 0 : setBoxDerivatives( virial[0] );
147 0 : }
148 :
149 10 : void DihedralCorrelation::calculateCV( const unsigned& mode, const std::vector<double>& masses, const std::vector<double>& charges,
150 : const std::vector<Vector>& pos, std::vector<double>& vals, std::vector<std::vector<Vector> >& derivs,
151 : std::vector<Tensor>& virial, const ActionAtomistic* aa ) {
152 10 : const Vector d10=delta(pos[1],pos[0]);
153 10 : const Vector d11=delta(pos[2],pos[1]);
154 10 : const Vector d12=delta(pos[3],pos[2]);
155 :
156 10 : Vector dd10,dd11,dd12;
157 : PLMD::Torsion t1;
158 10 : const double phi1 = t1.compute(d10,d11,d12,dd10,dd11,dd12);
159 :
160 10 : const Vector d20=delta(pos[5],pos[4]);
161 10 : const Vector d21=delta(pos[6],pos[5]);
162 10 : const Vector d22=delta(pos[7],pos[6]);
163 :
164 10 : Vector dd20,dd21,dd22;
165 : PLMD::Torsion t2;
166 10 : const double phi2 = t2.compute( d20, d21, d22, dd20, dd21, dd22 );
167 :
168 : // Calculate value
169 10 : const double diff = phi2 - phi1;
170 10 : vals[0] = 0.5*(1.+std::cos(diff));
171 : // Derivatives wrt phi1
172 10 : const double dval = 0.5*std::sin(diff);
173 10 : dd10 *= dval;
174 10 : dd11 *= dval;
175 10 : dd12 *= dval;
176 : // And add
177 10 : derivs[0][0]=dd10;
178 10 : derivs[0][1]=dd11-dd10;
179 10 : derivs[0][2]=dd12-dd11;
180 10 : derivs[0][3]=-dd12;
181 : // Derivative wrt phi2
182 10 : dd20 *= -dval;
183 10 : dd21 *= -dval;
184 10 : dd22 *= -dval;
185 : // And add
186 10 : derivs[0][4]=dd20;
187 10 : derivs[0][5]=dd21-dd20;
188 10 : derivs[0][6]=dd22-dd21;
189 10 : derivs[0][7]=-dd22;
190 10 : virial[0] = -(extProduct(d10,dd10)+extProduct(d11,dd11)+extProduct(d12,dd12)) - (extProduct(d20,dd20)+extProduct(d21,dd21)+extProduct(d22,dd22));
191 10 : }
192 :
193 : }
194 : }
|
