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 : #ifndef __PLUMED_colvar_Plane_h
23 : #define __PLUMED_colvar_Plane_h
24 : #include "Colvar.h"
25 : #include "ColvarInput.h"
26 : #include "tools/Pbc.h"
27 :
28 : #include <string>
29 : #include <cmath>
30 :
31 : namespace PLMD {
32 : namespace colvar {
33 :
34 : template <typename T>
35 : class Plane : 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 Plane(const ActionOptions&);
44 : static void parseAtomList( const int& num, std::vector<AtomNumber>& t, ActionAtomistic* aa );
45 : static unsigned getModeAndSetupValues( ActionWithValue* av );
46 : // active methods:
47 : void calculate() override;
48 : static void calculateCV( const ColvarInput<T>& cvin, ColvarOutput<T>& cvout );
49 : };
50 :
51 : template <typename T>
52 31 : void Plane<T>::registerKeywords( Keywords& keys ) {
53 31 : Colvar::registerKeywords( keys );
54 62 : keys.setDisplayName("PLANE");
55 31 : keys.add("atoms","ATOMS","the three or four atoms whose plane we are computing");
56 62 : keys.addOutputComponent("x","default","scalar/vector","the x-component of the vector that is normal to the plane containing the atoms");
57 62 : keys.addOutputComponent("y","default","scalar/vector","the y-component of the vector that is normal to the plane containing the atoms");
58 62 : keys.addOutputComponent("z","default","scalar/vector","the z-component of the vector that is normal to the plane containing the atoms");
59 31 : keys.add("hidden","NO_ACTION_LOG","suppresses printing from action on the log");
60 62 : keys.reset_style("NUMERICAL_DERIVATIVES","hidden");
61 31 : }
62 :
63 : template <typename T>
64 17 : void Plane<T>::parseAtomList( const int& num, std::vector<AtomNumber>& atoms, ActionAtomistic* aa ) {
65 34 : aa->parseAtomList("ATOMS",num,atoms);
66 17 : if(atoms.size()==3) {
67 14 : aa->log.printf(" containing atoms %d %d %d\n",atoms[0].serial(),atoms[1].serial(),atoms[2].serial());
68 14 : atoms.resize(4);
69 14 : atoms[3]=atoms[2];
70 14 : atoms[2]=atoms[1];
71 3 : } else if(atoms.size()==4) {
72 0 : aa->log.printf(" containing lines %d-%d and %d-%d\n",atoms[0].serial(),atoms[1].serial(),atoms[2].serial(),atoms[3].serial());
73 3 : } else if( num<0 || atoms.size()>0 ) {
74 0 : aa->error("Number of specified atoms should be either 3 or 4");
75 : }
76 17 : }
77 :
78 : template <typename T>
79 5 : unsigned Plane<T>::getModeAndSetupValues( ActionWithValue* av ) {
80 10 : av->addComponentWithDerivatives("x");
81 5 : av->componentIsNotPeriodic("x");
82 10 : av->addComponentWithDerivatives("y");
83 5 : av->componentIsNotPeriodic("y");
84 10 : av->addComponentWithDerivatives("z");
85 5 : av->componentIsNotPeriodic("z");
86 5 : return 0;
87 : }
88 :
89 : template <typename T>
90 2 : Plane<T>::Plane(const ActionOptions&ao):
91 : PLUMED_COLVAR_INIT(ao),
92 2 : pbc(true),
93 2 : value(3) {
94 2 : std::vector<AtomNumber> atoms;
95 2 : parseAtomList(-1,atoms,this);
96 2 : bool nopbc=!pbc;
97 2 : parseFlag("NOPBC",nopbc);
98 2 : pbc=!nopbc;
99 :
100 2 : if(pbc) {
101 2 : log.printf(" using periodic boundary conditions\n");
102 : } else {
103 0 : log.printf(" without periodic boundary conditions\n");
104 : }
105 :
106 2 : /*auto mode =*/ getModeAndSetupValues( this );
107 2 : requestAtoms(atoms);
108 2 : checkRead();
109 2 : }
110 :
111 : template <typename T>
112 115 : void Plane<T>::calculate() {
113 :
114 115 : if(pbc) {
115 115 : makeWhole();
116 : }
117 115 : auto cvout = ColvarOutput<double>::createColvarOutput(value,derivs,this);
118 115 : Plane<double>::calculateCV( ColvarInput<double>::createColvarInput( 0, getPositions(), this ), cvout );
119 115 : Value* valuex=getPntrToComponent("x");
120 115 : Value* valuey=getPntrToComponent("y");
121 115 : Value* valuez=getPntrToComponent("z");
122 :
123 575 : for(unsigned i=0; i<getPositions().size(); ++i) {
124 460 : setAtomsDerivatives( valuex, i, cvout.getAtomDerivatives(0,i) );
125 : }
126 115 : setBoxDerivatives( valuex, cvout.virial[0] );
127 115 : valuex->set( value[0] );
128 :
129 575 : for(unsigned i=0; i<getPositions().size(); ++i) {
130 460 : setAtomsDerivatives( valuey, i, cvout.getAtomDerivatives(1,i) );
131 : }
132 115 : setBoxDerivatives( valuey, cvout.virial[1] );
133 115 : valuey->set( value[1] );
134 :
135 575 : for(unsigned i=0; i<getPositions().size(); ++i) {
136 460 : setAtomsDerivatives( valuez, i, cvout.getAtomDerivatives(2,i) );
137 : }
138 115 : setBoxDerivatives( valuez, cvout.virial[2] );
139 115 : valuez->set( value[2] );
140 115 : }
141 :
142 : template <typename T>
143 301 : void Plane<T>::calculateCV( const ColvarInput<T>& cvin, ColvarOutput<T>& cvout ) {
144 301 : const auto d1=delta( cvin.pos[1], cvin.pos[0] );
145 301 : const auto d2=delta( cvin.pos[2], cvin.pos[3] );
146 301 : const auto cp = crossProduct( d1, d2 );
147 : using T33 = TensorTyped<T,3,3>;
148 301 : cvout.derivs[0][0] = crossProduct( Versors::xm<T>, d2 );
149 301 : cvout.derivs[0][1] = crossProduct( Versors::xp<T>, d2 );
150 301 : cvout.derivs[0][2] = crossProduct( Versors::xm<T>, d1 );
151 301 : cvout.derivs[0][3] = crossProduct( Versors::xp<T>, d1 );
152 602 : cvout.virial.set( 0, T33(d1,crossProduct(Versors::xp<T>, d2))
153 301 : + T33(d2, crossProduct(Versors::xm<T>, d1)));
154 301 : cvout.values[0] = cp[0];
155 :
156 301 : cvout.derivs[1][0] = crossProduct( Versors::ym<T>, d2 );
157 301 : cvout.derivs[1][1] = crossProduct( Versors::yp<T>, d2 );
158 301 : cvout.derivs[1][2] = crossProduct( Versors::ym<T>, d1 );
159 301 : cvout.derivs[1][3] = crossProduct( Versors::yp<T>, d1 );
160 602 : cvout.virial.set(1, T33(d1,crossProduct(Versors::yp<T>, d2))
161 301 : + T33(d2, crossProduct(Versors::ym<T>, d1)));
162 301 : cvout.values[1] = cp[1];
163 :
164 301 : cvout.derivs[2][0] = crossProduct( Versors::zm<T>, d2 );
165 301 : cvout.derivs[2][1] = crossProduct( Versors::zp<T>, d2 );
166 301 : cvout.derivs[2][2] = crossProduct( Versors::zm<T>, d1 );
167 301 : cvout.derivs[2][3] = crossProduct( Versors::zp<T>, d1 );
168 602 : cvout.virial.set(2, T33(d1,crossProduct(Versors::zp<T>, d2))
169 301 : + T33(d2, crossProduct(Versors::zm<T>, d1)));
170 301 : cvout.values[2] = cp[2];
171 301 : }
172 :
173 : }
174 : }
175 :
176 : #endif // __PLUMED_colvar_Plane_h
|
