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 : /// Get box shape
116 : const Tensor & getBox()const;
117 : /// Get the array of all positions
118 : const std::vector<Vector> & getPositions()const;
119 : /// Get the array of all masses
120 : const std::vector<double>& getMasses()const;
121 : /// Get the array of all charges
122 : const std::vector<double>& getCharges( const bool allowempty=false )const;
123 : /// Get the virial that is acting
124 : Tensor getVirial() const ;
125 : /// Get energy
126 : const double & getEnergy()const;
127 : /// Get mass of i-th atom
128 : double getMass(int i)const;
129 : /// Get charge of i-th atom
130 : double getCharge(int i)const;
131 : /// Get the force acting on a particular atom
132 : Vector getForce( const std::pair<std::size_t, std::size_t>& a ) const ;
133 : /// Add force to an atom
134 : void addForce( const std::pair<std::size_t, std::size_t>& a, const Vector& f );
135 : /// Get a reference to force on energy
136 : double & modifyForceOnEnergy();
137 : /// Get number of available atoms
138 : unsigned getNumberOfAtoms()const {
139 32693060 : return indexes.size();
140 : }
141 : /// Compute the pbc distance between two positions
142 : Vector pbcDistance(const Vector&,const Vector&)const;
143 : /// Applies PBCs to a seriens of positions or distances
144 : void pbcApply(std::vector<Vector>& dlist, unsigned max_index=0) const;
145 : /// Get the vector of absolute indexes
146 : virtual const std::vector<AtomNumber> & getAbsoluteIndexes()const;
147 : /// Get the absolute index of an atom
148 : AtomNumber getAbsoluteIndex(int i)const;
149 : /// Parse a list of atoms without a numbered keyword
150 : void parseAtomList(const std::string&key,std::vector<AtomNumber> &t);
151 : /// Parse an list of atom with a numbred keyword
152 : void parseAtomList(const std::string&key,const int num, std::vector<AtomNumber> &t);
153 : /// Interpret the atom selection. Just a wrapper to the static function with four arguments called interpretAtomList that passes xpos and this.
154 : void interpretAtomList( std::vector<std::string>& strings, std::vector<AtomNumber> &t);
155 : /// Convert a set of read in strings into an atom list (this is used in parseAtomList)
156 : static void interpretAtomList( std::vector<std::string>& strings, const std::vector<Value*>& xpos, Action* action, std::vector<AtomNumber> &t);
157 : /// This gets std::vector that contain the PLMD::Value objects that contain xpositions, ypositions, zpositions, masses and charges
158 : 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 );
159 : /// Change the box shape
160 : void changeBox( const Tensor& newbox );
161 : /// Get reference to Pbc
162 : const Pbc & getPbc() const;
163 : /// Add the forces to the atoms
164 : void setForcesOnAtoms( const std::vector<double>& forcesToApply, unsigned& ind );
165 : /// Add the virial forces
166 : void setForcesOnCell(const std::vector<double>& forcesToApply, unsigned& ind);
167 : /// Add the virial forces (span-like syntax)
168 : void setForcesOnCell(const double* forcesToApply, std::size_t size, unsigned& ind);
169 : /// Skip atom retrieval - use with care.
170 : /// If this function is called during initialization, then atoms are
171 : /// not going to be retrieved. Can be used for optimization. Notice that
172 : /// calling getPosition(int) in an Action where DoNotRetrieve() was called might
173 : /// lead to undefined behavior.
174 : void doNotRetrieve() {
175 82 : donotretrieve=true;
176 : }
177 : /// Skip atom forces - use with care.
178 : /// If this function is called during initialization, then forces are
179 : /// not going to be propagated. Can be used for optimization.
180 : void doNotForce() {
181 82 : donotforce=true;
182 : }
183 : /// Make atoms whole, assuming they are in the proper order
184 : void makeWhole();
185 : public:
186 :
187 : // virtual functions:
188 :
189 : explicit ActionAtomistic(const ActionOptions&ao);
190 : ~ActionAtomistic();
191 : static void registerKeywords( Keywords& keys );
192 :
193 : /// N.B. only pass an ActionWithValue to this routine if you know exactly what you
194 : /// are doing. The default will be correct for the vast majority of cases
195 : void calculateNumericalDerivatives( ActionWithValue* a=NULL ) override;
196 : /// Numerical derivative routine to use when using Actions that inherit from BOTH
197 : /// ActionWithArguments and ActionAtomistic
198 : void calculateAtomicNumericalDerivatives( ActionWithValue* a, const unsigned& startnum );
199 :
200 : virtual void retrieveAtoms( const bool& force=false );
201 : void lockRequests() override;
202 : void unlockRequests() override;
203 : const std::vector<AtomNumber> & getUnique()const;
204 : const std::vector<AtomNumber> & getUniqueLocal()const;
205 : /// Read in an input file containing atom positions and calculate the action for the atomic
206 : /// configuration therin
207 : void readAtomsFromPDB( const PDB& pdb ) override;
208 : /// Transfer the gradients
209 : void getGradient( const unsigned& ind, Vector& deriv, std::map<AtomNumber,Vector>& gradients ) const ;
210 473755 : ActionAtomistic* castToActionAtomistic() noexcept final {
211 473755 : return this;
212 : }
213 : virtual bool actionHasForces();
214 : };
215 :
216 : inline
217 : const Vector & ActionAtomistic::getPosition(int i)const {
218 127023996 : return actionPositions[i];
219 : }
220 :
221 : inline
222 3631283 : double ActionAtomistic::getMass(int i)const {
223 3631283 : if( !massesWereSet ) {
224 126432 : log.printf("WARNING: masses were not passed to plumed\n");
225 : }
226 3631283 : return masses[i];
227 : }
228 :
229 : inline
230 2905489 : double ActionAtomistic::getCharge(int i) const {
231 2905489 : if( !chargesWereSet ) {
232 0 : error("charges were not passed to plumed");
233 : }
234 2905489 : return charges[i];
235 : }
236 :
237 : inline
238 157 : const std::vector<AtomNumber> & ActionAtomistic::getAbsoluteIndexes()const {
239 157 : return indexes;
240 : }
241 :
242 : inline
243 : AtomNumber ActionAtomistic::getAbsoluteIndex(int i)const {
244 10009510 : return indexes[i];
245 : }
246 :
247 : inline
248 : const std::vector<Vector> & ActionAtomistic::getPositions()const {
249 492837 : return actionPositions;
250 : }
251 :
252 : inline
253 : const std::vector<double> & ActionAtomistic::getMasses()const {
254 : return masses;
255 : }
256 :
257 : inline
258 121311 : const std::vector<double> & ActionAtomistic::getCharges( const bool allowempty )const {
259 121311 : if( !allowempty && !chargesWereSet ) {
260 0 : error("charges were not passed to plumed");
261 : }
262 121311 : return charges;
263 : }
264 :
265 : inline
266 : const double & ActionAtomistic::getEnergy()const {
267 : return energy;
268 : }
269 :
270 : inline
271 : const Tensor & ActionAtomistic::getBox()const {
272 6449 : return actionPbc.getBox();
273 : }
274 :
275 : inline
276 : double & ActionAtomistic::modifyForceOnEnergy() {
277 : return forceOnEnergy;
278 : }
279 :
280 : inline
281 : const Pbc & ActionAtomistic::getPbc() const {
282 189438 : return actionPbc;
283 : }
284 :
285 : inline
286 207946 : void ActionAtomistic::lockRequests() {
287 450276 : lockRequestAtoms=true;
288 207946 : }
289 :
290 : inline
291 207946 : void ActionAtomistic::unlockRequests() {
292 450276 : lockRequestAtoms=false;
293 207946 : }
294 :
295 : inline
296 : const std::vector<AtomNumber> & ActionAtomistic::getUnique()const {
297 11456 : return unique;
298 : }
299 :
300 : inline
301 : const std::vector<AtomNumber> & ActionAtomistic::getUniqueLocal() const {
302 107066 : return unique_local;
303 : }
304 :
305 : inline
306 928875 : Vector ActionAtomistic::getGlobalPosition(const std::pair<std::size_t,std::size_t>& a) const {
307 : Vector pos;
308 928875 : pos[0]=xpos[a.first]->data[a.second];
309 928875 : pos[1]=ypos[a.first]->data[a.second];
310 928875 : pos[2]=zpos[a.first]->data[a.second];
311 928875 : return pos;
312 : }
313 :
314 : inline
315 909912 : void ActionAtomistic::setGlobalPosition(const std::pair<std::size_t, std::size_t>& a, const Vector& pos ) {
316 909912 : xpos[a.first]->data[a.second]=pos[0];
317 909912 : ypos[a.first]->data[a.second]=pos[1];
318 909912 : zpos[a.first]->data[a.second]=pos[2];
319 909912 : }
320 :
321 : inline
322 15886 : Vector ActionAtomistic::getForce( const std::pair<std::size_t, std::size_t>& a ) const {
323 : Vector f;
324 15886 : f[0]=xpos[a.first]->getForce(a.second);
325 15886 : f[1]=ypos[a.first]->getForce(a.second);
326 15886 : f[2]=zpos[a.first]->getForce(a.second);
327 15886 : return f;
328 : }
329 :
330 : inline
331 9946 : void ActionAtomistic::addForce( const std::pair<std::size_t, std::size_t>& a, const Vector& f ) {
332 9946 : xpos[a.first]->addForce( a.second, f[0] );
333 9946 : ypos[a.first]->addForce( a.second, f[1] );
334 9946 : zpos[a.first]->addForce( a.second, f[2] );
335 9946 : }
336 :
337 : inline
338 : Vector ActionAtomistic::pbcDistance(const Vector &v1,const Vector &v2)const {
339 29741892 : return actionPbc.distance(v1,v2);
340 : }
341 :
342 : }
343 : #endif
|
