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