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 "core/ActionRegister.h"
25 : #include "MultiColvarTemplate.h"
26 : #include "tools/Pbc.h"
27 :
28 : #include <string>
29 : #include <cmath>
30 :
31 : //+PLUMEDOC COLVAR PLANE
32 : /*
33 : Calculate the plane perpendicular to two vectors in order to represent the orientation of a planar molecule.
34 :
35 : To calculate the orientation of the plane connecting atoms 1, 2 and 3 you use an input like this:
36 :
37 : ```plumed
38 : p: PLANE ATOMS=1,2,3
39 : PRINT ARG=p.x,p.y,p.z FILE=colvar
40 : ```
41 :
42 : The three components, p.x, p.y and p.z, output by the PLANE action here are the x, y and z components of the normal
43 : vector to the plane that is obtained by taking the cross product between the vector connecting atoms 1 and 2 and
44 : the vector connecting atoms 2 and 3.
45 :
46 : To calculate the cross product of the vector connecting atoms 1 and 2 and the vector connecting atoms 3 and 4 you use
47 : an input like this:
48 :
49 : ```plumed
50 : p: PLANE ATOMS=1,2,3,4
51 : PRINT ARG=p.x,p.y,p.z FILE=colvar
52 : ```
53 :
54 : If you have multiple molecules and wish to determine the orientations of the planes containing all them with one line of PLUMED
55 : input you can use an input like this:
56 :
57 : ```plumed
58 : p: PLANE ATOMS1=1,2,3 ATOMS2=4,5,6 ATOMS3=7,8,9 ATOMS4=10,11,12
59 : PRINT ARG=p.x,p.y,p.z FILE=colvar
60 : ```
61 :
62 : The output from this command consists of 3 vectors with 4 components. These vectors, p.x, p.y and p.z, contain the x, y and z components
63 : of the normals to the planes of the molecules. Commands similar to this are useful for variables that can be used to monitor
64 : nucleation of molecular crystals such as [SMAC](SMAC.md).
65 :
66 : */
67 : //+ENDPLUMEDOC
68 :
69 : //+PLUMEDOC COLVAR PLANE_SCALAR
70 : /*
71 : Calculate the plane perpendicular to two vectors in order to represent the orientation of a planar molecule.
72 :
73 : \par Examples
74 :
75 : */
76 : //+ENDPLUMEDOC
77 :
78 : //+PLUMEDOC COLVAR PLANE_VECTOR
79 : /*
80 : Calculate the plane perpendicular to two vectors in order to represent the orientation of a planar molecule multiple times.
81 :
82 : \par Examples
83 :
84 : */
85 : //+ENDPLUMEDOC
86 :
87 : namespace PLMD {
88 : namespace colvar {
89 :
90 : class Plane : public Colvar {
91 : private:
92 : bool pbc;
93 : std::vector<double> value, masses, charges;
94 : std::vector<std::vector<Vector> > derivs;
95 : std::vector<Tensor> virial;
96 : public:
97 : static void registerKeywords( Keywords& keys );
98 : explicit Plane(const ActionOptions&);
99 : static void parseAtomList( const int& num, std::vector<AtomNumber>& t, ActionAtomistic* aa );
100 : static unsigned getModeAndSetupValues( ActionWithValue* av );
101 : // active methods:
102 : void calculate() override;
103 : static void calculateCV( const unsigned& mode, const std::vector<double>& masses, const std::vector<double>& charges,
104 : const std::vector<Vector>& pos, std::vector<double>& vals, std::vector<std::vector<Vector> >& derivs,
105 : std::vector<Tensor>& virial, const ActionAtomistic* aa );
106 : };
107 :
108 : typedef ColvarShortcut<Plane> PlaneShortcut;
109 : PLUMED_REGISTER_ACTION(PlaneShortcut,"PLANE")
110 : PLUMED_REGISTER_ACTION(Plane,"PLANE_SCALAR")
111 : typedef MultiColvarTemplate<Plane> PlaneMulti;
112 : PLUMED_REGISTER_ACTION(PlaneMulti,"PLANE_VECTOR")
113 :
114 10 : void Plane::registerKeywords( Keywords& keys ) {
115 10 : Colvar::registerKeywords( keys );
116 20 : keys.setDisplayName("PLANE");
117 10 : keys.add("atoms","ATOMS","the three or four atoms whose plane we are computing");
118 20 : keys.addOutputComponent("x","default","scalar/vector","the x-component of the vector that is normal to the plane containing the atoms");
119 20 : keys.addOutputComponent("y","default","scalar/vector","the y-component of the vector that is normal to the plane containing the atoms");
120 20 : keys.addOutputComponent("z","default","scalar/vector","the z-component of the vector that is normal to the plane containing the atoms");
121 10 : keys.add("hidden","NO_ACTION_LOG","suppresses printing from action on the log");
122 10 : }
123 :
124 9 : void Plane::parseAtomList( const int& num, std::vector<AtomNumber>& atoms, ActionAtomistic* aa ) {
125 18 : aa->parseAtomList("ATOMS",num,atoms);
126 9 : if(atoms.size()==3) {
127 8 : aa->log.printf(" containing atoms %d %d %d\n",atoms[0].serial(),atoms[1].serial(),atoms[2].serial());
128 8 : atoms.resize(4);
129 8 : atoms[3]=atoms[2];
130 8 : atoms[2]=atoms[1];
131 1 : } else if(atoms.size()==4) {
132 0 : aa->log.printf(" containing lines %d-%d and %d-%d\n",atoms[0].serial(),atoms[1].serial(),atoms[2].serial(),atoms[3].serial());
133 1 : } else if( num<0 || atoms.size()>0 ) {
134 0 : aa->error("Number of specified atoms should be either 3 or 4");
135 : }
136 9 : }
137 :
138 1 : unsigned Plane::getModeAndSetupValues( ActionWithValue* av ) {
139 2 : av->addComponentWithDerivatives("x");
140 1 : av->componentIsNotPeriodic("x");
141 2 : av->addComponentWithDerivatives("y");
142 1 : av->componentIsNotPeriodic("y");
143 2 : av->addComponentWithDerivatives("z");
144 1 : av->componentIsNotPeriodic("z");
145 1 : return 0;
146 : }
147 :
148 0 : Plane::Plane(const ActionOptions&ao):
149 : PLUMED_COLVAR_INIT(ao),
150 0 : pbc(true),
151 0 : value(3),
152 0 : derivs(3),
153 0 : virial(3) {
154 0 : for(unsigned i=0; i<3; ++i) {
155 0 : derivs[i].resize(4);
156 : }
157 : std::vector<AtomNumber> atoms;
158 0 : parseAtomList(-1,atoms,this);
159 0 : bool nopbc=!pbc;
160 0 : parseFlag("NOPBC",nopbc);
161 0 : pbc=!nopbc;
162 :
163 0 : if(pbc) {
164 0 : log.printf(" using periodic boundary conditions\n");
165 : } else {
166 0 : log.printf(" without periodic boundary conditions\n");
167 : }
168 :
169 0 : unsigned mode = getModeAndSetupValues( this );
170 0 : requestAtoms(atoms);
171 0 : checkRead();
172 0 : }
173 :
174 0 : void Plane::calculate() {
175 :
176 0 : if(pbc) {
177 0 : makeWhole();
178 : }
179 0 : calculateCV( 0, masses, charges, getPositions(), value, derivs, virial, this );
180 0 : setValue( value[0] );
181 0 : for(unsigned i=0; i<derivs[0].size(); ++i) {
182 0 : setAtomsDerivatives( i, derivs[0][i] );
183 : }
184 0 : setBoxDerivatives( virial[0] );
185 0 : }
186 :
187 176 : void Plane::calculateCV( const unsigned& mode, const std::vector<double>& masses, const std::vector<double>& charges,
188 : const std::vector<Vector>& pos, std::vector<double>& vals, std::vector<std::vector<Vector> >& derivs,
189 : std::vector<Tensor>& virial, const ActionAtomistic* aa ) {
190 176 : Vector d1=delta( pos[1], pos[0] );
191 176 : Vector d2=delta( pos[2], pos[3] );
192 176 : Vector cp = crossProduct( d1, d2 );
193 :
194 176 : derivs[0][0] = crossProduct( Vector(-1.0,0,0), d2 );
195 176 : derivs[0][1] = crossProduct( Vector(+1.0,0,0), d2 );
196 176 : derivs[0][2] = crossProduct( Vector(-1.0,0,0), d1 );
197 176 : derivs[0][3] = crossProduct( Vector(+1.0,0,0), d1 );
198 176 : virial[0] = Tensor(d1,crossProduct(Vector(+1.0,0,0), d2)) + Tensor( d2, crossProduct(Vector(-1.0,0,0), d1));
199 176 : vals[0] = cp[0];
200 :
201 176 : derivs[1][0] = crossProduct( Vector(0,-1.0,0), d2 );
202 176 : derivs[1][1] = crossProduct( Vector(0,+1.0,0), d2 );
203 176 : derivs[1][2] = crossProduct( Vector(0,-1.0,0), d1 );
204 176 : derivs[1][3] = crossProduct( Vector(0,+1.0,0), d1 );
205 176 : virial[1] = Tensor(d1,crossProduct(Vector(0,+1.0,0), d2)) + Tensor( d2, crossProduct(Vector(0,-1.0,0), d1));
206 176 : vals[1] = cp[1];
207 :
208 176 : derivs[2][0] = crossProduct( Vector(0,0,-1.0), d2 );
209 176 : derivs[2][1] = crossProduct( Vector(0,0,+1.0), d2 );
210 176 : derivs[2][2] = crossProduct( Vector(0,0,-1.0), d1 );
211 176 : derivs[2][3] = crossProduct( Vector(0,0,+1.0), d1 );
212 176 : virial[2] = Tensor(d1,crossProduct(Vector(0,0,+1.0), d2)) + Tensor( d2, crossProduct(Vector(0,0,-1.0), d1));
213 176 : vals[2] = cp[2];
214 176 : }
215 :
216 : }
217 : }
218 :
219 :
220 :
|
