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 : #ifndef __PLUMED_colvar_DihedralCorrelation_h
23 : #define __PLUMED_colvar_DihedralCorrelation_h
24 : #include "Colvar.h"
25 : #include "ColvarInput.h"
26 : #include "tools/Torsion.h"
27 :
28 : #include <string>
29 : #include <cmath>
30 :
31 : namespace PLMD {
32 : namespace colvar {
33 :
34 : template <typename T>
35 : class DihedralCorrelation : public Colvar {
36 : private:
37 : bool pbc;
38 : std::vector<double> value;
39 : std::vector<double> derivs;
40 : public:
41 : using precision=T;
42 : static void registerKeywords( Keywords& keys );
43 : explicit DihedralCorrelation(const ActionOptions&);
44 : static void parseAtomList( const int& num, std::vector<AtomNumber>& t, ActionAtomistic* aa );
45 : static unsigned getModeAndSetupValues( ActionWithValue* av );
46 : void calculate() override;
47 : static void calculateCV( const ColvarInput<T>& cvin, ColvarOutput<T>& cvout );
48 : };
49 :
50 : template <typename T>
51 19 : void DihedralCorrelation<T>::registerKeywords( Keywords& keys ) {
52 19 : Colvar::registerKeywords( keys );
53 38 : keys.setDisplayName("DIHEDRAL_CORRELATION");
54 19 : keys.add("atoms","ATOMS","the set of 8 atoms that are being used to calculate this quantity");
55 19 : keys.add("hidden","NO_ACTION_LOG","suppresses printing from action on the log");
56 38 : keys.setValueDescription("scalar/vector","the DIHEDRAL_CORRELATION for these atoms");
57 38 : keys.reset_style("NUMERICAL_DERIVATIVES","hidden");
58 19 : }
59 :
60 : template <typename T>
61 1 : DihedralCorrelation<T>::DihedralCorrelation(const ActionOptions&ao):
62 : PLUMED_COLVAR_INIT(ao),
63 1 : pbc(true),
64 1 : value(1) {
65 1 : std::vector<AtomNumber> atoms;
66 1 : parseAtomList(-1,atoms,this);
67 1 : if( atoms.size()!=8 ) {
68 0 : error("Number of specified atoms should be 8");
69 : }
70 :
71 1 : bool nopbc=!pbc;
72 1 : parseFlag("NOPBC",nopbc);
73 1 : pbc=!nopbc;
74 :
75 1 : if(pbc) {
76 1 : log.printf(" using periodic boundary conditions\n");
77 : } else {
78 0 : log.printf(" without periodic boundary conditions\n");
79 : }
80 1 : addValueWithDerivatives();
81 1 : setNotPeriodic();
82 1 : requestAtoms(atoms);
83 1 : checkRead();
84 1 : }
85 :
86 : template <typename T>
87 4 : void DihedralCorrelation<T>::parseAtomList( const int& num, std::vector<AtomNumber>& t, ActionAtomistic* aa ) {
88 4 : aa->parseAtomList("ATOMS",num,t);
89 4 : if( num<0 && t.size()!=8 ) {
90 0 : aa->error("Number of specified atoms should be 8");
91 : }
92 4 : if( t.size()==8 ) {
93 3 : aa->log.printf(" correlation between dihedral angle for atoms %d %d %d %d and atoms %d %d %d %d\n",
94 : t[0].serial(),t[1].serial(),t[2].serial(),t[3].serial(),t[4].serial(),t[5].serial(),t[6].serial(),t[7].serial());
95 : }
96 4 : }
97 :
98 : template <typename T>
99 1 : unsigned DihedralCorrelation<T>::getModeAndSetupValues( ActionWithValue* av ) {
100 1 : av->addValueWithDerivatives();
101 1 : av->setNotPeriodic();
102 1 : return 0;
103 : }
104 :
105 : template <typename T>
106 4 : void DihedralCorrelation<T>::calculate() {
107 :
108 4 : if(pbc) {
109 4 : makeWhole();
110 : }
111 4 : auto cvout = ColvarOutput<double>::createColvarOutput(value,derivs,this);
112 4 : DihedralCorrelation<double>::calculateCV( ColvarInput<double>::createColvarInput( 0, getPositions(), this ), cvout );
113 4 : setValue( value[0] );
114 36 : for(unsigned i=0; i<getPositions().size(); ++i) {
115 32 : setAtomsDerivatives( i, cvout.getAtomDerivatives(0, i) );
116 : }
117 4 : setBoxDerivatives( cvout.virial[0] );
118 4 : }
119 :
120 : template <typename T>
121 24 : void DihedralCorrelation<T>::calculateCV( const ColvarInput<T>& cvin, ColvarOutput<T>& cvout ) {
122 24 : const auto d10=delta(cvin.pos[1],cvin.pos[0]);
123 24 : const auto d11=delta(cvin.pos[2],cvin.pos[1]);
124 24 : const auto d12=delta(cvin.pos[3],cvin.pos[2]);
125 :
126 24 : VectorTyped<T,3> dd10,dd11,dd12;
127 24 : const T phi1 = PLMD::Torsion::compute(d10,d11,d12,dd10,dd11,dd12);
128 :
129 24 : const auto d20=delta(cvin.pos[5],cvin.pos[4]);
130 24 : const auto d21=delta(cvin.pos[6],cvin.pos[5]);
131 24 : const auto d22=delta(cvin.pos[7],cvin.pos[6]);
132 :
133 24 : VectorTyped<T,3> dd20,dd21,dd22;
134 24 : const T phi2 = PLMD::Torsion::compute( d20, d21, d22, dd20, dd21, dd22 );
135 :
136 : // Calculate value
137 24 : const T diff = phi2 - phi1;
138 24 : cvout.values[0] = 0.5*(1.+std::cos(diff));
139 : // Derivatives wrt phi1
140 24 : const T dval = 0.5*std::sin(diff);
141 : dd10 *= dval;
142 : dd11 *= dval;
143 : dd12 *= dval;
144 : // And add
145 : cvout.derivs[0][0]=dd10;
146 : cvout.derivs[0][1]=dd11-dd10;
147 : cvout.derivs[0][2]=dd12-dd11;
148 : cvout.derivs[0][3]=-dd12;
149 : // Derivative wrt phi2
150 24 : dd20 *= -dval;
151 : dd21 *= -dval;
152 : dd22 *= -dval;
153 : // And add
154 : cvout.derivs[0][4]=dd20;
155 : cvout.derivs[0][5]=dd21-dd20;
156 : cvout.derivs[0][6]=dd22-dd21;
157 : cvout.derivs[0][7]=-dd22;
158 48 : cvout.virial.set( 0,
159 24 : -(extProduct(d10,dd10)+extProduct(d11,dd11)+extProduct(d12,dd12))
160 24 : -(extProduct(d20,dd20)+extProduct(d21,dd21)+extProduct(d22,dd22)) );
161 24 : }
162 :
163 : } // namespace colvar
164 : } // namespace PLMD
165 : #endif //__PLUMED_colvar_DihedralCorrelation_h
|
