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 10 : void Plane::registerKeywords( Keywords& keys ) {
87 10 : Colvar::registerKeywords( keys ); keys.setDisplayName("PLANE");
88 20 : keys.add("atoms","ATOMS","the three or four atoms whose plane we are computing");
89 20 : keys.addOutputComponent("x","default","scalar/vector","the x-component of the vector that is normal to the plane containing the atoms");
90 20 : keys.addOutputComponent("y","default","scalar/vector","the y-component of the vector that is normal to the plane containing the atoms");
91 20 : keys.addOutputComponent("z","default","scalar/vector","the z-component of the vector that is normal to the plane containing the atoms");
92 20 : keys.add("hidden","NO_ACTION_LOG","suppresses printing from action on the log");
93 10 : }
94 :
95 9 : void Plane::parseAtomList( const int& num, std::vector<AtomNumber>& atoms, ActionAtomistic* aa ) {
96 18 : aa->parseAtomList("ATOMS",num,atoms);
97 9 : if(atoms.size()==3) {
98 8 : aa->log.printf(" containing atoms %d %d %d\n",atoms[0].serial(),atoms[1].serial(),atoms[2].serial());
99 8 : atoms.resize(4);
100 8 : atoms[3]=atoms[2];
101 8 : atoms[2]=atoms[1];
102 1 : } else if(atoms.size()==4) {
103 0 : aa->log.printf(" containing lines %d-%d and %d-%d\n",atoms[0].serial(),atoms[1].serial(),atoms[2].serial(),atoms[3].serial());
104 1 : } else if( num<0 || atoms.size()>0 ) aa->error("Number of specified atoms should be either 3 or 4");
105 9 : }
106 :
107 1 : unsigned Plane::getModeAndSetupValues( ActionWithValue* av ) {
108 2 : av->addComponentWithDerivatives("x"); av->componentIsNotPeriodic("x");
109 2 : av->addComponentWithDerivatives("y"); av->componentIsNotPeriodic("y");
110 2 : av->addComponentWithDerivatives("z"); av->componentIsNotPeriodic("z");
111 1 : return 0;
112 : }
113 :
114 0 : Plane::Plane(const ActionOptions&ao):
115 : PLUMED_COLVAR_INIT(ao),
116 0 : pbc(true),
117 0 : value(3),
118 0 : derivs(3),
119 0 : virial(3)
120 : {
121 0 : for(unsigned i=0; i<3; ++i) derivs[i].resize(4);
122 0 : std::vector<AtomNumber> atoms; parseAtomList(-1,atoms,this);
123 0 : bool nopbc=!pbc;
124 0 : parseFlag("NOPBC",nopbc);
125 0 : pbc=!nopbc;
126 :
127 0 : if(pbc) log.printf(" using periodic boundary conditions\n");
128 0 : else log.printf(" without periodic boundary conditions\n");
129 :
130 0 : unsigned mode = getModeAndSetupValues( this );
131 0 : requestAtoms(atoms);
132 0 : checkRead();
133 0 : }
134 :
135 0 : void Plane::calculate() {
136 :
137 0 : if(pbc) makeWhole();
138 0 : calculateCV( 0, masses, charges, getPositions(), value, derivs, virial, this );
139 0 : setValue( value[0] );
140 0 : for(unsigned i=0; i<derivs[0].size(); ++i) setAtomsDerivatives( i, derivs[0][i] );
141 0 : setBoxDerivatives( virial[0] );
142 0 : }
143 :
144 176 : void Plane::calculateCV( const unsigned& mode, const std::vector<double>& masses, const std::vector<double>& charges,
145 : const std::vector<Vector>& pos, std::vector<double>& vals, std::vector<std::vector<Vector> >& derivs,
146 : std::vector<Tensor>& virial, const ActionAtomistic* aa ) {
147 176 : Vector d1=delta( pos[1], pos[0] );
148 176 : Vector d2=delta( pos[2], pos[3] );
149 176 : Vector cp = crossProduct( d1, d2 );
150 :
151 176 : derivs[0][0] = crossProduct( Vector(-1.0,0,0), d2 );
152 176 : derivs[0][1] = crossProduct( Vector(+1.0,0,0), d2 );
153 176 : derivs[0][2] = crossProduct( Vector(-1.0,0,0), d1 );
154 176 : derivs[0][3] = crossProduct( Vector(+1.0,0,0), d1 );
155 176 : virial[0] = Tensor(d1,crossProduct(Vector(+1.0,0,0), d2)) + Tensor( d2, crossProduct(Vector(-1.0,0,0), d1));
156 176 : vals[0] = cp[0];
157 :
158 176 : derivs[1][0] = crossProduct( Vector(0,-1.0,0), d2 );
159 176 : derivs[1][1] = crossProduct( Vector(0,+1.0,0), d2 );
160 176 : derivs[1][2] = crossProduct( Vector(0,-1.0,0), d1 );
161 176 : derivs[1][3] = crossProduct( Vector(0,+1.0,0), d1 );
162 176 : virial[1] = Tensor(d1,crossProduct(Vector(0,+1.0,0), d2)) + Tensor( d2, crossProduct(Vector(0,-1.0,0), d1));
163 176 : vals[1] = cp[1];
164 :
165 176 : derivs[2][0] = crossProduct( Vector(0,0,-1.0), d2 );
166 176 : derivs[2][1] = crossProduct( Vector(0,0,+1.0), d2 );
167 176 : derivs[2][2] = crossProduct( Vector(0,0,-1.0), d1 );
168 176 : derivs[2][3] = crossProduct( Vector(0,0,+1.0), d1 );
169 176 : virial[2] = Tensor(d1,crossProduct(Vector(0,0,+1.0), d2)) + Tensor( d2, crossProduct(Vector(0,0,-1.0), d1));
170 176 : vals[2] = cp[2];
171 176 : }
172 :
173 : }
174 : }
175 :
176 :
177 :
|
