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