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 ); keys.setDisplayName("DIHEDRAL_CORRELATION");
89 20 : keys.add("atoms","ATOMS","the set of 8 atoms that are being used to calculate this quantity");
90 20 : keys.add("hidden","NO_ACTION_LOG","suppresses printing from action on the log");
91 10 : keys.setValueDescription("the DIHEDRAL_CORRELATION for these atoms");
92 10 : }
93 :
94 0 : DihedralCorrelation::DihedralCorrelation(const ActionOptions&ao):
95 : PLUMED_COLVAR_INIT(ao),
96 0 : pbc(true),
97 0 : value(1),
98 0 : derivs(1),
99 0 : virial(1)
100 : {
101 0 : derivs[0].resize(8); std::vector<AtomNumber> atoms;
102 0 : parseAtomList(-1,atoms,this);
103 0 : if( atoms.size()!=8 ) error("Number of specified atoms should be 8");
104 :
105 0 : bool nopbc=!pbc;
106 0 : parseFlag("NOPBC",nopbc);
107 0 : pbc=!nopbc;
108 :
109 0 : if(pbc) log.printf(" using periodic boundary conditions\n");
110 0 : else log.printf(" without periodic boundary conditions\n");
111 0 : }
112 :
113 3 : void DihedralCorrelation::parseAtomList( const int& num, std::vector<AtomNumber>& t, ActionAtomistic* aa ) {
114 3 : aa->parseAtomList("ATOMS",num,t);
115 3 : if( num<0 && t.size()!=8 ) aa->error("Number of specified atoms should be 8");
116 3 : if( t.size()==8 ) {
117 2 : aa->log.printf(" correlation between dihedral angle for atoms %d %d %d %d and atoms %d %d %d %d\n",
118 : t[0].serial(),t[1].serial(),t[2].serial(),t[3].serial(),t[4].serial(),t[5].serial(),t[6].serial(),t[7].serial());
119 : }
120 3 : }
121 :
122 1 : unsigned DihedralCorrelation::getModeAndSetupValues( ActionWithValue* av ) {
123 2 : av->addValueWithDerivatives(); av->setNotPeriodic(); return 0;
124 : }
125 :
126 0 : void DihedralCorrelation::calculate() {
127 :
128 0 : if(pbc) makeWhole();
129 0 : calculateCV( 0, masses, charges, getPositions(), value, derivs, virial, this );
130 0 : setValue( value[0] );
131 0 : for(unsigned i=0; i<derivs[0].size(); ++i) setAtomsDerivatives( i, derivs[0][i] );
132 0 : setBoxDerivatives( virial[0] );
133 0 : }
134 :
135 10 : void DihedralCorrelation::calculateCV( const unsigned& mode, const std::vector<double>& masses, const std::vector<double>& charges,
136 : const std::vector<Vector>& pos, std::vector<double>& vals, std::vector<std::vector<Vector> >& derivs,
137 : std::vector<Tensor>& virial, const ActionAtomistic* aa ) {
138 10 : const Vector d10=delta(pos[1],pos[0]);
139 10 : const Vector d11=delta(pos[2],pos[1]);
140 10 : const Vector d12=delta(pos[3],pos[2]);
141 :
142 10 : Vector dd10,dd11,dd12;
143 : PLMD::Torsion t1;
144 10 : const double phi1 = t1.compute(d10,d11,d12,dd10,dd11,dd12);
145 :
146 10 : const Vector d20=delta(pos[5],pos[4]);
147 10 : const Vector d21=delta(pos[6],pos[5]);
148 10 : const Vector d22=delta(pos[7],pos[6]);
149 :
150 10 : Vector dd20,dd21,dd22;
151 : PLMD::Torsion t2;
152 10 : const double phi2 = t2.compute( d20, d21, d22, dd20, dd21, dd22 );
153 :
154 : // Calculate value
155 10 : const double diff = phi2 - phi1;
156 10 : vals[0] = 0.5*(1.+std::cos(diff));
157 : // Derivatives wrt phi1
158 10 : const double dval = 0.5*std::sin(diff);
159 10 : dd10 *= dval;
160 10 : dd11 *= dval;
161 10 : dd12 *= dval;
162 : // And add
163 10 : derivs[0][0]=dd10;
164 10 : derivs[0][1]=dd11-dd10;
165 10 : derivs[0][2]=dd12-dd11;
166 10 : derivs[0][3]=-dd12;
167 : // Derivative wrt phi2
168 10 : dd20 *= -dval;
169 10 : dd21 *= -dval;
170 10 : dd22 *= -dval;
171 : // And add
172 10 : derivs[0][4]=dd20;
173 10 : derivs[0][5]=dd21-dd20;
174 10 : derivs[0][6]=dd22-dd21;
175 10 : derivs[0][7]=-dd22;
176 10 : virial[0] = -(extProduct(d10,dd10)+extProduct(d11,dd11)+extProduct(d12,dd12)) - (extProduct(d20,dd20)+extProduct(d21,dd21)+extProduct(d22,dd22));
177 10 : }
178 :
179 : }
180 : }
|
