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 : \par Examples
36 :
37 :
38 : */
39 : //+ENDPLUMEDOC
40 :
41 : //+PLUMEDOC COLVAR PLANE_SCALAR
42 : /*
43 : Calculate the plane perpendicular to two vectors in order to represent the orientation of a planar molecule.
44 :
45 : \par Examples
46 :
47 : */
48 : //+ENDPLUMEDOC
49 :
50 : //+PLUMEDOC COLVAR PLANE_VECTOR
51 : /*
52 : Calculate the plane perpendicular to two vectors in order to represent the orientation of a planar molecule multiple times.
53 :
54 : \par Examples
55 :
56 : */
57 : //+ENDPLUMEDOC
58 :
59 : namespace PLMD {
60 : namespace colvar {
61 :
62 : class Plane : public Colvar {
63 : private:
64 : bool pbc;
65 : std::vector<double> value, masses, charges;
66 : std::vector<std::vector<Vector> > derivs;
67 : std::vector<Tensor> virial;
68 : public:
69 : static void registerKeywords( Keywords& keys );
70 : explicit Plane(const ActionOptions&);
71 : static void parseAtomList( const int& num, std::vector<AtomNumber>& t, ActionAtomistic* aa );
72 : static unsigned getModeAndSetupValues( ActionWithValue* av );
73 : // active methods:
74 : void calculate() override;
75 : static void calculateCV( const unsigned& mode, const std::vector<double>& masses, const std::vector<double>& charges,
76 : const std::vector<Vector>& pos, std::vector<double>& vals, std::vector<std::vector<Vector> >& derivs,
77 : std::vector<Tensor>& virial, const ActionAtomistic* aa );
78 : };
79 :
80 : typedef ColvarShortcut<Plane> PlaneShortcut;
81 : PLUMED_REGISTER_ACTION(PlaneShortcut,"PLANE")
82 : PLUMED_REGISTER_ACTION(Plane,"PLANE_SCALAR")
83 : typedef MultiColvarTemplate<Plane> PlaneMulti;
84 : PLUMED_REGISTER_ACTION(PlaneMulti,"PLANE_VECTOR")
85 :
86 18 : void Plane::registerKeywords( Keywords& keys ) {
87 18 : Colvar::registerKeywords( keys );
88 18 : keys.setDisplayName("PLANE");
89 36 : keys.add("atoms","ATOMS","the three or four atoms whose plane we are computing");
90 36 : keys.addOutputComponent("x","default","the x-component of the vector that is normal to the plane containing the atoms");
91 36 : keys.addOutputComponent("y","default","the y-component of the vector that is normal to the plane containing the atoms");
92 36 : keys.addOutputComponent("z","default","the z-component of the vector that is normal to the plane containing the atoms");
93 36 : keys.add("hidden","NO_ACTION_LOG","suppresses printing from action on the log");
94 18 : }
95 :
96 11 : void Plane::parseAtomList( const int& num, std::vector<AtomNumber>& atoms, ActionAtomistic* aa ) {
97 22 : aa->parseAtomList("ATOMS",num,atoms);
98 11 : if(atoms.size()==3) {
99 10 : aa->log.printf(" containing atoms %d %d %d\n",atoms[0].serial(),atoms[1].serial(),atoms[2].serial());
100 10 : atoms.resize(4);
101 10 : atoms[3]=atoms[2];
102 10 : atoms[2]=atoms[1];
103 1 : } else if(atoms.size()==4) {
104 0 : aa->log.printf(" containing lines %d-%d and %d-%d\n",atoms[0].serial(),atoms[1].serial(),atoms[2].serial(),atoms[3].serial());
105 1 : } else if( num<0 || atoms.size()>0 ) {
106 0 : aa->error("Number of specified atoms should be either 3 or 4");
107 : }
108 11 : }
109 :
110 3 : unsigned Plane::getModeAndSetupValues( ActionWithValue* av ) {
111 6 : av->addComponentWithDerivatives("x");
112 3 : av->componentIsNotPeriodic("x");
113 6 : av->addComponentWithDerivatives("y");
114 3 : av->componentIsNotPeriodic("y");
115 6 : av->addComponentWithDerivatives("z");
116 3 : av->componentIsNotPeriodic("z");
117 3 : return 0;
118 : }
119 :
120 2 : Plane::Plane(const ActionOptions&ao):
121 : PLUMED_COLVAR_INIT(ao),
122 2 : pbc(true),
123 2 : value(3),
124 2 : derivs(3),
125 4 : virial(3) {
126 8 : for(unsigned i=0; i<3; ++i) {
127 6 : derivs[i].resize(4);
128 : }
129 : std::vector<AtomNumber> atoms;
130 2 : parseAtomList(-1,atoms,this);
131 2 : bool nopbc=!pbc;
132 2 : parseFlag("NOPBC",nopbc);
133 2 : pbc=!nopbc;
134 :
135 2 : if(pbc) {
136 2 : log.printf(" using periodic boundary conditions\n");
137 : } else {
138 0 : log.printf(" without periodic boundary conditions\n");
139 : }
140 :
141 2 : unsigned mode = getModeAndSetupValues( this );
142 2 : requestAtoms(atoms);
143 2 : checkRead();
144 2 : }
145 :
146 115 : void Plane::calculate() {
147 :
148 115 : if(pbc) {
149 115 : makeWhole();
150 : }
151 115 : calculateCV( 0, masses, charges, getPositions(), value, derivs, virial, this );
152 115 : Value* valuex=getPntrToComponent("x");
153 115 : Value* valuey=getPntrToComponent("y");
154 115 : Value* valuez=getPntrToComponent("z");
155 :
156 575 : for(unsigned i=0; i<derivs[0].size(); ++i) {
157 460 : setAtomsDerivatives( valuex, i, derivs[0][i] );
158 : }
159 115 : setBoxDerivatives( valuex, virial[0] );
160 115 : valuex->set( value[0] );
161 :
162 575 : for(unsigned i=0; i<derivs[1].size(); ++i) {
163 460 : setAtomsDerivatives( valuey, i, derivs[1][i] );
164 : }
165 115 : setBoxDerivatives( valuey, virial[1] );
166 115 : valuey->set( value[1] );
167 :
168 575 : for(unsigned i=0; i<derivs[2].size(); ++i) {
169 460 : setAtomsDerivatives( valuez, i, derivs[2][i] );
170 : }
171 115 : setBoxDerivatives( valuez, virial[2] );
172 115 : valuez->set( value[2] );
173 115 : }
174 :
175 291 : void Plane::calculateCV( const unsigned& mode, const std::vector<double>& masses, const std::vector<double>& charges,
176 : const std::vector<Vector>& pos, std::vector<double>& vals, std::vector<std::vector<Vector> >& derivs,
177 : std::vector<Tensor>& virial, const ActionAtomistic* aa ) {
178 291 : Vector d1=delta( pos[1], pos[0] );
179 291 : Vector d2=delta( pos[2], pos[3] );
180 291 : Vector cp = crossProduct( d1, d2 );
181 :
182 291 : derivs[0][0] = crossProduct( Vector(-1.0,0,0), d2 );
183 291 : derivs[0][1] = crossProduct( Vector(+1.0,0,0), d2 );
184 291 : derivs[0][2] = crossProduct( Vector(-1.0,0,0), d1 );
185 291 : derivs[0][3] = crossProduct( Vector(+1.0,0,0), d1 );
186 291 : virial[0] = Tensor(d1,crossProduct(Vector(+1.0,0,0), d2)) + Tensor( d2, crossProduct(Vector(-1.0,0,0), d1));
187 291 : vals[0] = cp[0];
188 :
189 291 : derivs[1][0] = crossProduct( Vector(0,-1.0,0), d2 );
190 291 : derivs[1][1] = crossProduct( Vector(0,+1.0,0), d2 );
191 291 : derivs[1][2] = crossProduct( Vector(0,-1.0,0), d1 );
192 291 : derivs[1][3] = crossProduct( Vector(0,+1.0,0), d1 );
193 291 : virial[1] = Tensor(d1,crossProduct(Vector(0,+1.0,0), d2)) + Tensor( d2, crossProduct(Vector(0,-1.0,0), d1));
194 291 : vals[1] = cp[1];
195 :
196 291 : derivs[2][0] = crossProduct( Vector(0,0,-1.0), d2 );
197 291 : derivs[2][1] = crossProduct( Vector(0,0,+1.0), d2 );
198 291 : derivs[2][2] = crossProduct( Vector(0,0,-1.0), d1 );
199 291 : derivs[2][3] = crossProduct( Vector(0,0,+1.0), d1 );
200 291 : virial[2] = Tensor(d1,crossProduct(Vector(0,0,+1.0), d2)) + Tensor( d2, crossProduct(Vector(0,0,-1.0), d1));
201 291 : vals[2] = cp[2];
202 291 : }
203 :
204 : }
205 : }
206 :
207 :
208 :
|
