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