Line data Source code
1 : /* +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
2 : Copyright (c) 2011-2023 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 "core/ActionRegister.h"
26 : #include "tools/Angle.h"
27 :
28 : namespace PLMD {
29 : namespace colvar {
30 :
31 : //+PLUMEDOC COLVAR ANGLE
32 : /*
33 : Calculate an angle.
34 :
35 : This command can be used to compute the angle between three atoms. Alternatively
36 : if four atoms appear in the atom
37 : specification it calculates the angle between
38 : two vectors identified by two pairs of atoms.
39 :
40 : If _three_ atoms are given, the angle is defined as:
41 : \f[
42 : \theta=\arccos\left(\frac{ {\bf r}_{21}\cdot {\bf r}_{23}}{
43 : |{\bf r}_{21}| |{\bf r}_{23}|}\right)
44 : \f]
45 : Here \f$ {\bf r}_{ij}\f$ is the distance vector among the
46 : \f$i\f$th and the \f$j\f$th listed atom.
47 :
48 : If _four_ atoms are given, the angle is defined as:
49 : \f[
50 : \theta=\arccos\left(\frac{ {\bf r}_{21}\cdot {\bf r}_{34}}{
51 : |{\bf r}_{21}| |{\bf r}_{34}|}\right)
52 : \f]
53 :
54 : Notice that angles defined in this way are non-periodic variables and
55 : their value is limited by definition between 0 and \f$\pi\f$.
56 :
57 : The vectors \f$ {\bf r}_{ij}\f$ are by default evaluated taking
58 : periodic boundary conditions into account.
59 : This behavior can be changed with the NOPBC flag.
60 :
61 : \par Examples
62 :
63 : This command tells plumed to calculate the angle between the vector connecting atom 1 to atom 2 and
64 : the vector connecting atom 2 to atom 3 and to print it on file COLVAR1. At the same time,
65 : the angle between vector connecting atom 1 to atom 2 and the vector connecting atom 3 to atom 4 is printed
66 : on file COLVAR2.
67 : \plumedfile
68 :
69 : a: ANGLE ATOMS=1,2,3
70 : # equivalently one could state:
71 : # a: ANGLE ATOMS=1,2,2,3
72 :
73 : b: ANGLE ATOMS=1,2,3,4
74 :
75 : PRINT ARG=a FILE=COLVAR1
76 : PRINT ARG=b FILE=COLVAR2
77 : \endplumedfile
78 :
79 :
80 : */
81 : //+ENDPLUMEDOC
82 :
83 : //+PLUMEDOC COLVAR ANGLE_SCALAR
84 : /*
85 : Calculate an angle.
86 :
87 : \par Examples
88 :
89 : */
90 : //+ENDPLUMEDOC
91 :
92 : //+PLUMEDOC COLVAR ANGLE_VECTOR
93 : /*
94 : Calculate multiple angles.
95 :
96 : \par Examples
97 :
98 : */
99 : //+ENDPLUMEDOC
100 :
101 : class Angle : public Colvar {
102 : bool pbc;
103 : std::vector<double> value, masses, charges;
104 : std::vector<std::vector<Vector> > derivs;
105 : std::vector<Tensor> virial;
106 : public:
107 : explicit Angle(const ActionOptions&);
108 : // active methods:
109 : void calculate() override;
110 : static void registerKeywords( Keywords& keys );
111 : static void parseAtomList( const int& num, std::vector<AtomNumber>& t, ActionAtomistic* aa );
112 : static unsigned getModeAndSetupValues( ActionWithValue* av );
113 : static void calculateCV( const unsigned& mode, const std::vector<double>& masses, const std::vector<double>& charges,
114 : const std::vector<Vector>& pos, std::vector<double>& vals, std::vector<std::vector<Vector> >& derivs,
115 : std::vector<Tensor>& virial, const ActionAtomistic* aa );
116 : };
117 :
118 : typedef ColvarShortcut<Angle> AngleShortcut;
119 : PLUMED_REGISTER_ACTION(AngleShortcut,"ANGLE")
120 : PLUMED_REGISTER_ACTION(Angle,"ANGLE_SCALAR")
121 : typedef MultiColvarTemplate<Angle> AngleMulti;
122 : PLUMED_REGISTER_ACTION(AngleMulti,"ANGLE_VECTOR")
123 :
124 89 : void Angle::registerKeywords( Keywords& keys ) {
125 89 : Colvar::registerKeywords(keys); keys.setDisplayName("ANGLE");
126 178 : keys.add("atoms","ATOMS","the list of atoms involved in this collective variable (either 3 or 4 atoms)");
127 178 : keys.add("hidden","NO_ACTION_LOG","suppresses printing from action on the log");
128 89 : keys.setValueDescription("the ANGLE involving these atoms");
129 89 : }
130 :
131 73 : void Angle::parseAtomList( const int& num, std::vector<AtomNumber>& atoms, ActionAtomistic* aa ) {
132 146 : aa->parseAtomList("ATOMS",num,atoms);
133 73 : if(atoms.size()==3) {
134 62 : aa->log.printf(" between atoms %d %d %d\n",atoms[0].serial(),atoms[1].serial(),atoms[2].serial());
135 62 : atoms.resize(4);
136 62 : atoms[3]=atoms[2];
137 62 : atoms[2]=atoms[1];
138 11 : } else if(atoms.size()==4) {
139 7 : aa->log.printf(" between lines %d-%d and %d-%d\n",atoms[0].serial(),atoms[1].serial(),atoms[2].serial(),atoms[3].serial());
140 4 : } else if( num<0 || atoms.size()>0 ) aa->error("Number of specified atoms should be either 3 or 4");
141 72 : }
142 :
143 3 : unsigned Angle::getModeAndSetupValues( ActionWithValue* av ) {
144 6 : av->addValueWithDerivatives(); av->setNotPeriodic(); return 0;
145 : }
146 :
147 23 : Angle::Angle(const ActionOptions&ao):
148 : PLUMED_COLVAR_INIT(ao),
149 23 : pbc(true),
150 23 : value(1),
151 24 : derivs(1),
152 46 : virial(1)
153 : {
154 23 : derivs[0].resize(4);
155 23 : std::vector<AtomNumber> atoms; parseAtomList( -1, atoms, this );
156 22 : bool nopbc=!pbc;
157 22 : parseFlag("NOPBC",nopbc);
158 22 : pbc=!nopbc;
159 :
160 22 : if(pbc) log.printf(" using periodic boundary conditions\n");
161 0 : else log.printf(" without periodic boundary conditions\n");
162 :
163 45 : addValueWithDerivatives(); setNotPeriodic();
164 22 : requestAtoms(atoms);
165 22 : checkRead();
166 25 : }
167 :
168 : // calculator
169 319 : void Angle::calculate() {
170 :
171 319 : if(pbc) makeWhole();
172 319 : calculateCV( 0, masses, charges, getPositions(), value, derivs, virial, this );
173 319 : setValue( value[0] );
174 1595 : for(unsigned i=0; i<derivs[0].size(); ++i) setAtomsDerivatives( i, derivs[0][i] );
175 319 : setBoxDerivatives( virial[0] );
176 319 : }
177 :
178 554 : void Angle::calculateCV( const unsigned& mode, const std::vector<double>& masses, const std::vector<double>& charges,
179 : const std::vector<Vector>& pos, std::vector<double>& vals, std::vector<std::vector<Vector> >& derivs,
180 : std::vector<Tensor>& virial, const ActionAtomistic* aa ) {
181 554 : Vector dij,dik;
182 554 : dij=delta(pos[2],pos[3]);
183 554 : dik=delta(pos[1],pos[0]);
184 554 : Vector ddij,ddik; PLMD::Angle a;
185 554 : vals[0]=a.compute(dij,dik,ddij,ddik);
186 554 : derivs[0][0]=ddik; derivs[0][1]=-ddik;
187 554 : derivs[0][2]=-ddij; derivs[0][3]=ddij;
188 554 : setBoxDerivativesNoPbc( pos, derivs, virial );
189 554 : }
190 :
191 : }
192 : }
193 :
194 :
195 :
|
