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_core_ActionAtomistic_h
23 : #define __PLUMED_core_ActionAtomistic_h
24 :
25 : #include "Action.h"
26 : #include "tools/Tensor.h"
27 : #include "tools/Pbc.h"
28 : #include "tools/ForwardDecl.h"
29 : #include "Value.h"
30 : #include <vector>
31 : #include <map>
32 :
33 : namespace PLMD {
34 :
35 : class Pbc;
36 : class PDB;
37 : class GenericMolInfo;
38 : class Tree;
39 :
40 : namespace colvar {
41 : class SelectMassCharge;
42 : }
43 :
44 : /// \ingroup MULTIINHERIT
45 : /// Action used to create objects that access the positions of the atoms from the MD code
46 : class ActionAtomistic :
47 : virtual public Action {
48 : friend class Group;
49 : friend class DomainDecomposition;
50 : friend class colvar::SelectMassCharge;
51 : friend class ActionWithVirtualAtom;
52 :
53 : std::vector<AtomNumber> indexes; // the set of needed atoms
54 : std::vector<std::size_t> value_depends; // The list of values that are being used
55 : std::vector<std::pair<std::size_t, std::size_t > > atom_value_ind; // The list of values and indices for the atoms that are being used
56 : std::vector<std::pair<std::size_t,std::vector<std::size_t>>> atom_value_ind_grouped;
57 : /// unique should be an ordered set since we later create a vector containing the corresponding indexes
58 : std::vector<AtomNumber> unique;
59 : /// unique_local should be an ordered set since we later create a vector containing the corresponding indexes
60 : bool unique_local_needs_update;
61 : std::vector<AtomNumber> unique_local;
62 : std::vector<Vector> actionPositions; // positions of the needed atoms
63 : double energy;
64 : Value* boxValue;
65 : ForwardDecl<Pbc> pbc_fwd;
66 : Pbc& actionPbc=*pbc_fwd;
67 : std::vector<double> masses;
68 : std::vector<double> charges;
69 :
70 : std::vector<Vector> forces; // forces on the needed atoms
71 : double forceOnEnergy;
72 :
73 : double forceOnExtraCV;
74 :
75 : bool lockRequestAtoms; // forbid changes to request atoms
76 :
77 : bool donotretrieve;
78 : bool donotforce;
79 :
80 : // EMST
81 : GenericMolInfo* actionMoldat{nullptr};
82 : std::unique_ptr<Tree> tree;
83 :
84 : /// Values that hold information about atom positions and charges
85 : std::vector<Value*> xpos, ypos, zpos, masv, chargev;
86 : void updateUniqueLocal( const bool& useunique, const std::vector<int>& g2l );
87 : protected:
88 : bool massesWereSet;
89 : bool chargesWereSet;
90 : void setExtraCV(const std::string &name);
91 : /// Used to interpret whether this index is a virtual atom or a real atom
92 : std::pair<std::size_t, std::size_t> getValueIndices( const AtomNumber& i ) const ;
93 : public:
94 : /// Request an array of atoms.
95 : /// This method is used to ask for a list of atoms. Atoms
96 : /// should be asked for by number. If this routine is called
97 : /// during the simulation, atoms will be available at the next step
98 : /// MAYBE WE HAVE TO FIND SOMETHING MORE CLEAR FOR DYNAMIC
99 : /// LISTS OF ATOMS
100 : void requestAtoms(const std::vector<AtomNumber> & a, const bool clearDep=true);
101 : /// Get position of i-th atom (access by relative index)
102 : const Vector & getPosition(int)const;
103 : /// Get position of i-th atom (access by absolute AtomNumber).
104 : /// With direct access to the global atom array.
105 : /// \warning Should be only used by actions that need to read the shared position array.
106 : /// This array is insensitive to local changes such as makeWhole(), numerical derivatives, etc.
107 : Vector getGlobalPosition(const std::pair<std::size_t,std::size_t>& ) const ;
108 : /// Modify position of i-th atom (access by absolute AtomNumber).
109 : /// \warning Should be only used by actions that need to modify the shared position array.
110 : /// This array is insensitive to local changes such as makeWhole(), numerical derivatives, etc.
111 : void setGlobalPosition(const std::pair<std::size_t,std::size_t>&, const Vector& pos);
112 : /// Get total number of atoms, including virtual ones.
113 : /// Can be used to make a loop on modifyGlobalPosition or getGlobalPosition.
114 : unsigned getTotAtoms()const;
115 : /// Translate all the atoms, including virtual ones, in the shared position arrays.
116 : /// These functions loop directly on the arrays, which is much faster than calling
117 : /// getGlobalPosition and setGlobalPosition for every atom.
118 : /// \warning Should be only used by actions that need to modify the shared position array.
119 : void translateAllGlobalPositions(const Vector& shift);
120 : /// Rotate all the atoms, including virtual ones, in the shared position arrays: pos -> matmul(rotation,pos)
121 : void rotateAllGlobalPositions(const Tensor& rotation);
122 : /// Move all the atoms, including virtual ones, in the shared position arrays: pos -> matmul(rotation,pos-shiftBefore)+shiftAfter
123 : void rotateAllGlobalPositions(const Tensor& rotation, const Vector& shiftBefore, const Vector& shiftAfter);
124 : /// Get the sum of the forces on all the atoms, including virtual ones
125 : Vector getTotalForceOnAllAtoms() const;
126 : /// Rotate the forces on all the atoms, including virtual ones: f -> matmul(rotation,f).
127 : /// Returns the sum of the rotated forces.
128 : Vector rotateForcesOnAllAtoms(const Tensor& rotation);
129 : /// Get box shape
130 : const Tensor & getBox()const;
131 : /// Get the array of all positions
132 : const std::vector<Vector> & getPositions()const;
133 : /// Get the array of all masses
134 : const std::vector<double>& getMasses()const;
135 : /// Get the array of all charges
136 : const std::vector<double>& getCharges( const bool allowempty=false )const;
137 : /// Get the virial that is acting
138 : Tensor getVirial() const ;
139 : /// Get energy
140 : const double & getEnergy()const;
141 : /// Get mass of i-th atom
142 : double getMass(int i)const;
143 : /// Get charge of i-th atom
144 : double getCharge(int i)const;
145 : /// Get the force acting on a particular atom
146 : Vector getForce( const std::pair<std::size_t, std::size_t>& a ) const ;
147 : /// Add force to an atom
148 : void addForce( const std::pair<std::size_t, std::size_t>& a, const Vector& f );
149 : /// Get a reference to force on energy
150 : double & modifyForceOnEnergy();
151 : /// Get number of available atoms
152 : unsigned getNumberOfAtoms()const {
153 32426149 : return indexes.size();
154 : }
155 : /// Compute the pbc distance between two positions
156 : Vector pbcDistance(const Vector&,const Vector&)const;
157 : /// Applies PBCs to a seriens of positions or distances
158 : void pbcApply(std::vector<Vector>& dlist, unsigned max_index=0) const;
159 : /// Get the vector of absolute indexes
160 : virtual const std::vector<AtomNumber> & getAbsoluteIndexes()const;
161 : /// Get the absolute index of an atom
162 : AtomNumber getAbsoluteIndex(int i)const;
163 : /// Parse a list of atoms without a numbered keyword
164 : void parseAtomList(const std::string&key,std::vector<AtomNumber> &t);
165 : /// Parse an list of atom with a numbred keyword
166 : void parseAtomList(const std::string&key,const int num, std::vector<AtomNumber> &t);
167 : /// Interpret the atom selection. Just a wrapper to the static function with four arguments called interpretAtomList that passes xpos and this.
168 : void interpretAtomList( std::vector<std::string>& strings, std::vector<AtomNumber> &t);
169 : /// Convert a set of read in strings into an atom list (this is used in parseAtomList)
170 : static void interpretAtomList( std::vector<std::string>& strings, const std::vector<Value*>& xpos, Action* action, std::vector<AtomNumber> &t);
171 : /// This gets std::vector that contain the PLMD::Value objects that contain xpositions, ypositions, zpositions, masses and charges
172 : static void getAtomValuesFromPlumedObject( const PlumedMain& plumed, std::vector<Value*>& xpos, std::vector<Value*>& ypos, std::vector<Value*>& zpos, std::vector<Value*>& masv, std::vector<Value*>& chargev );
173 : /// Change the box shape
174 : void changeBox( const Tensor& newbox );
175 : /// Get reference to Pbc
176 : const Pbc & getPbc() const;
177 : /// Add the forces to the atoms
178 : void setForcesOnAtoms( const std::vector<double>& forcesToApply, unsigned& ind );
179 : /// Add the virial forces
180 : void setForcesOnCell(const std::vector<double>& forcesToApply, unsigned& ind);
181 : /// Add the virial forces (span-like syntax)
182 : void setForcesOnCell(const double* forcesToApply, std::size_t size, unsigned& ind);
183 : /// Skip atom retrieval - use with care.
184 : /// If this function is called during initialization, then atoms are
185 : /// not going to be retrieved. Can be used for optimization. Notice that
186 : /// calling getPosition(int) in an Action where DoNotRetrieve() was called might
187 : /// lead to undefined behavior.
188 : void doNotRetrieve() {
189 113 : donotretrieve=true;
190 : }
191 : /// Skip atom forces - use with care.
192 : /// If this function is called during initialization, then forces are
193 : /// not going to be propagated. Can be used for optimization.
194 : void doNotForce() {
195 113 : donotforce=true;
196 : }
197 : /// Make atoms whole, assuming they are in the proper order
198 : void makeWhole();
199 : public:
200 :
201 : // virtual functions:
202 :
203 : explicit ActionAtomistic(const ActionOptions&ao);
204 : ~ActionAtomistic();
205 : static void registerKeywords( Keywords& keys );
206 :
207 : /// N.B. only pass an ActionWithValue to this routine if you know exactly what you
208 : /// are doing. The default will be correct for the vast majority of cases
209 : void calculateNumericalDerivatives( ActionWithValue* a=NULL ) override;
210 : /// Numerical derivative routine to use when using Actions that inherit from BOTH
211 : /// ActionWithArguments and ActionAtomistic
212 : void calculateAtomicNumericalDerivatives( ActionWithValue* a, const unsigned& startnum );
213 :
214 : virtual void retrieveAtoms( const bool& force=false );
215 : void lockRequests() override;
216 : void unlockRequests() override;
217 : const std::vector<AtomNumber> & getUnique()const;
218 : const std::vector<AtomNumber> & getUniqueLocal()const;
219 : /// Read in an input file containing atom positions and calculate the action for the atomic
220 : /// configuration therin
221 : void readAtomsFromPDB( const PDB& pdb ) override;
222 : /// Transfer the gradients
223 : void getGradient( const unsigned& ind, Vector& deriv, std::map<AtomNumber,Vector>& gradients ) const ;
224 638994 : ActionAtomistic* castToActionAtomistic() noexcept final {
225 638994 : return this;
226 : }
227 : virtual bool actionHasForces();
228 : };
229 :
230 : inline
231 : const Vector & ActionAtomistic::getPosition(int i)const {
232 48642428 : return actionPositions[i];
233 : }
234 :
235 : inline
236 3641707 : double ActionAtomistic::getMass(int i)const {
237 3641707 : if( !massesWereSet ) {
238 126432 : log.printf("WARNING: masses were not passed to plumed\n");
239 : }
240 3641707 : return masses[i];
241 : }
242 :
243 : inline
244 3264061 : double ActionAtomistic::getCharge(int i) const {
245 3264061 : if( !chargesWereSet ) {
246 0 : error("charges were not passed to plumed");
247 : }
248 3264061 : return charges[i];
249 : }
250 :
251 : inline
252 163 : const std::vector<AtomNumber> & ActionAtomistic::getAbsoluteIndexes()const {
253 163 : return indexes;
254 : }
255 :
256 : inline
257 : AtomNumber ActionAtomistic::getAbsoluteIndex(int i)const {
258 9125263 : return indexes[i];
259 : }
260 :
261 : inline
262 : const std::vector<Vector> & ActionAtomistic::getPositions()const {
263 656390 : return actionPositions;
264 : }
265 :
266 : inline
267 : const std::vector<double> & ActionAtomistic::getMasses()const {
268 : return masses;
269 : }
270 :
271 : inline
272 294130 : const std::vector<double> & ActionAtomistic::getCharges( const bool allowempty )const {
273 294130 : if( !allowempty && !chargesWereSet ) {
274 0 : error("charges were not passed to plumed");
275 : }
276 294130 : return charges;
277 : }
278 :
279 : inline
280 : const double & ActionAtomistic::getEnergy()const {
281 : return energy;
282 : }
283 :
284 : inline
285 : const Tensor & ActionAtomistic::getBox()const {
286 6520 : return actionPbc.getBox();
287 : }
288 :
289 : inline
290 : double & ActionAtomistic::modifyForceOnEnergy() {
291 : return forceOnEnergy;
292 : }
293 :
294 : inline
295 : const Pbc & ActionAtomistic::getPbc() const {
296 434424 : return actionPbc;
297 : }
298 :
299 : inline
300 379403 : void ActionAtomistic::lockRequests() {
301 608696 : lockRequestAtoms=true;
302 379403 : }
303 :
304 : inline
305 379403 : void ActionAtomistic::unlockRequests() {
306 608696 : lockRequestAtoms=false;
307 379403 : }
308 :
309 : inline
310 : const std::vector<AtomNumber> & ActionAtomistic::getUnique()const {
311 6804 : return unique;
312 : }
313 :
314 : inline
315 : const std::vector<AtomNumber> & ActionAtomistic::getUniqueLocal() const {
316 16719 : return unique_local;
317 : }
318 :
319 : inline
320 1223200 : Vector ActionAtomistic::getGlobalPosition(const std::pair<std::size_t,std::size_t>& a) const {
321 : Vector pos;
322 1223200 : pos[0]=xpos[a.first]->data[a.second];
323 1223200 : pos[1]=ypos[a.first]->data[a.second];
324 1223200 : pos[2]=zpos[a.first]->data[a.second];
325 1223200 : return pos;
326 : }
327 :
328 : inline
329 1195441 : void ActionAtomistic::setGlobalPosition(const std::pair<std::size_t, std::size_t>& a, const Vector& pos ) {
330 1195441 : xpos[a.first]->data[a.second]=pos[0];
331 1195441 : ypos[a.first]->data[a.second]=pos[1];
332 1195441 : zpos[a.first]->data[a.second]=pos[2];
333 1195441 : }
334 :
335 : inline
336 : Vector ActionAtomistic::getForce( const std::pair<std::size_t, std::size_t>& a ) const {
337 : Vector f;
338 : f[0]=xpos[a.first]->getForce(a.second);
339 : f[1]=ypos[a.first]->getForce(a.second);
340 : f[2]=zpos[a.first]->getForce(a.second);
341 : return f;
342 : }
343 :
344 : inline
345 876 : void ActionAtomistic::addForce( const std::pair<std::size_t, std::size_t>& a, const Vector& f ) {
346 876 : xpos[a.first]->addForce( a.second, f[0] );
347 876 : ypos[a.first]->addForce( a.second, f[1] );
348 876 : zpos[a.first]->addForce( a.second, f[2] );
349 876 : }
350 :
351 : inline
352 : Vector ActionAtomistic::pbcDistance(const Vector &v1,const Vector &v2)const {
353 29506286 : return actionPbc.distance(v1,v2);
354 : }
355 :
356 : }
357 : #endif
|
