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 : #include "ActionAtomistic.h"
23 : #include "PlumedMain.h"
24 : #include "ActionSet.h"
25 : #include "GenericMolInfo.h"
26 : #include <vector>
27 : #include <string>
28 : #include "ActionWithValue.h"
29 : #include "Colvar.h"
30 : #include "ActionWithVirtualAtom.h"
31 : #include "tools/Exception.h"
32 : #include "Atoms.h"
33 : #include "tools/Pbc.h"
34 : #include "tools/PDB.h"
35 :
36 : namespace PLMD {
37 :
38 10380 : ActionAtomistic::~ActionAtomistic() {
39 : // forget the pending request
40 10380 : atoms.remove(this);
41 10380 : }
42 :
43 10380 : ActionAtomistic::ActionAtomistic(const ActionOptions&ao):
44 : Action(ao),
45 10380 : lockRequestAtoms(false),
46 10380 : donotretrieve(false),
47 10380 : donotforce(false),
48 10380 : atoms(plumed.getAtoms()) {
49 10380 : atoms.add(this);
50 : // if(atoms.getNatoms()==0) error("Cannot perform calculations involving atoms without atoms");
51 10380 : }
52 :
53 10982 : void ActionAtomistic::registerKeywords( Keywords& keys ) {
54 : (void) keys; // avoid warning
55 10982 : }
56 :
57 :
58 10371 : void ActionAtomistic::requestAtoms(const std::vector<AtomNumber> & a, const bool clearDep) {
59 10371 : plumed_massert(!lockRequestAtoms,"requested atom list can only be changed in the prepare() method");
60 10371 : int nat=a.size();
61 10371 : indexes=a;
62 10371 : positions.resize(nat);
63 10371 : forces.resize(nat);
64 10371 : masses.resize(nat);
65 10371 : charges.resize(nat);
66 10371 : int n=atoms.positions.size();
67 10371 : if(clearDep) {
68 10251 : clearDependencies();
69 : }
70 10371 : unique.clear();
71 1408480 : for(unsigned i=0; i<indexes.size(); i++) {
72 1398110 : if(indexes[i].index()>=n) {
73 : std::string num;
74 1 : Tools::convert( indexes[i].serial(),num );
75 3 : error("atom " + num + " out of range");
76 : }
77 1398109 : if(atoms.isVirtualAtom(indexes[i])) {
78 7225 : addDependency(atoms.getVirtualAtomsAction(indexes[i]));
79 : }
80 : // only real atoms are requested to lower level Atoms class
81 : else {
82 1390884 : unique.push_back(indexes[i]);
83 : }
84 : }
85 10370 : Tools::removeDuplicates(unique);
86 10370 : updateUniqueLocal();
87 10370 : atoms.unique.clear();
88 10370 : }
89 :
90 186514274 : Vector ActionAtomistic::pbcDistance(const Vector &v1,const Vector &v2)const {
91 186514274 : return pbc.distance(v1,v2);
92 : }
93 :
94 222049 : void ActionAtomistic::pbcApply(std::vector<Vector>& dlist, unsigned max_index)const {
95 222049 : pbc.apply(dlist, max_index);
96 222049 : }
97 :
98 195 : void ActionAtomistic::calculateNumericalDerivatives( ActionWithValue* a ) {
99 195 : calculateAtomicNumericalDerivatives( a, 0 );
100 195 : }
101 :
102 0 : void ActionAtomistic::changeBox( const Tensor& newbox ) {
103 0 : pbc.setBox( newbox );
104 0 : }
105 :
106 493 : void ActionAtomistic::calculateAtomicNumericalDerivatives( ActionWithValue* a, const unsigned& startnum ) {
107 493 : if(!a) {
108 280 : a=dynamic_cast<ActionWithValue*>(this);
109 280 : plumed_massert(a,"only Actions with a value can be differentiated");
110 : }
111 :
112 493 : const size_t nval=a->getNumberOfComponents();
113 493 : const size_t natoms=getNumberOfAtoms();
114 493 : std::vector<Vector> value(nval*natoms);
115 493 : std::vector<Tensor> valuebox(nval);
116 493 : std::vector<Vector> savedPositions(natoms);
117 : const double delta=std::sqrt(epsilon);
118 :
119 23936 : for(int i=0; i<natoms; i++)
120 93772 : for(int k=0; k<3; k++) {
121 70329 : savedPositions[i][k]=positions[i][k];
122 70329 : positions[i][k]=positions[i][k]+delta;
123 70329 : a->calculate();
124 70329 : positions[i][k]=savedPositions[i][k];
125 211947 : for(int j=0; j<nval; j++) {
126 141618 : value[j*natoms+i][k]=a->getOutputQuantity(j);
127 : }
128 : }
129 493 : Tensor box(pbc.getBox());
130 1972 : for(int i=0; i<3; i++)
131 5916 : for(int k=0; k<3; k++) {
132 4437 : double arg0=box(i,k);
133 215424 : for(int j=0; j<natoms; j++) {
134 210987 : positions[j]=pbc.realToScaled(positions[j]);
135 : }
136 4437 : box(i,k)=box(i,k)+delta;
137 4437 : pbc.setBox(box);
138 215424 : for(int j=0; j<natoms; j++) {
139 210987 : positions[j]=pbc.scaledToReal(positions[j]);
140 : }
141 4437 : a->calculate();
142 4437 : box(i,k)=arg0;
143 4437 : pbc.setBox(box);
144 215424 : for(int j=0; j<natoms; j++) {
145 210987 : positions[j]=savedPositions[j];
146 : }
147 19557 : for(int j=0; j<nval; j++) {
148 15120 : valuebox[j](i,k)=a->getOutputQuantity(j);
149 : }
150 : }
151 :
152 493 : a->calculate();
153 493 : a->clearDerivatives();
154 2173 : for(int j=0; j<nval; j++) {
155 1680 : Value* v=a->copyOutput(j);
156 : double ref=v->get();
157 1680 : if(v->hasDerivatives()) {
158 26188 : for(int i=0; i<natoms; i++)
159 101784 : for(int k=0; k<3; k++) {
160 76338 : double d=(value[j*natoms+i][k]-ref)/delta;
161 76338 : v->addDerivative(startnum+3*i+k,d);
162 : }
163 742 : Tensor virial;
164 2968 : for(int i=0; i<3; i++)
165 8904 : for(int k=0; k<3; k++) {
166 6678 : virial(i,k)= (valuebox[j](i,k)-ref)/delta;
167 : }
168 : // BE CAREFUL WITH NON ORTHOROMBIC CELL
169 742 : virial=-matmul(box.transpose(),virial);
170 2968 : for(int i=0; i<3; i++)
171 8904 : for(int k=0; k<3; k++) {
172 6678 : v->addDerivative(startnum+3*natoms+3*k+i,virial(k,i));
173 : }
174 : }
175 : }
176 493 : }
177 :
178 11905 : void ActionAtomistic::parseAtomList(const std::string&key, std::vector<AtomNumber> &t) {
179 11905 : parseAtomList(key,-1,t);
180 11904 : }
181 :
182 14415 : void ActionAtomistic::parseAtomList(const std::string&key,const int num, std::vector<AtomNumber> &t) {
183 14415 : plumed_massert( keywords.style(key,"atoms") || keywords.style(key,"hidden"), "keyword " + key + " should be registered as atoms");
184 : std::vector<std::string> strings;
185 14415 : if( num<0 ) {
186 11905 : parseVector(key,strings);
187 11904 : if(strings.empty()) {
188 : return;
189 : }
190 : } else {
191 2510 : if ( !parseNumberedVector(key,num,strings) ) {
192 : return;
193 : }
194 : }
195 11713 : interpretAtomList( strings, t );
196 14415 : }
197 :
198 11914 : void ActionAtomistic::interpretAtomList(std::vector<std::string>& strings, std::vector<AtomNumber> &t) {
199 11914 : Tools::interpretRanges(strings);
200 11914 : t.resize(0);
201 458942 : for(unsigned i=0; i<strings.size(); ++i) {
202 : AtomNumber atom;
203 447028 : bool ok=Tools::convertNoexcept(strings[i],atom); // this is converting strings to AtomNumbers
204 447028 : if(ok) {
205 438642 : t.push_back(atom);
206 : }
207 : // here we check if this is a special symbol for MOLINFO
208 447028 : if( !ok && strings[i].compare(0,1,"@")==0 ) {
209 763 : std::string symbol=strings[i].substr(1);
210 763 : if(symbol=="allatoms") {
211 10 : const auto n=plumed.getAtoms().getNatoms() + plumed.getAtoms().getNVirtualAtoms();
212 10 : t.reserve(t.size()+n);
213 365 : for(unsigned i=0; i<n; i++) {
214 355 : t.push_back(AtomNumber::index(i));
215 : }
216 : ok=true;
217 753 : } else if(symbol=="mdatoms") {
218 4 : const auto n=plumed.getAtoms().getNatoms();
219 4 : t.reserve(t.size()+n);
220 228 : for(unsigned i=0; i<n; i++) {
221 224 : t.push_back(AtomNumber::index(i));
222 : }
223 : ok=true;
224 : } else {
225 749 : auto* moldat=plumed.getActionSet().selectLatest<GenericMolInfo*>(this);
226 749 : if( moldat ) {
227 : std::vector<AtomNumber> atom_list;
228 749 : moldat->interpretSymbol( symbol, atom_list );
229 749 : if( atom_list.size()>0 ) {
230 : ok=true;
231 749 : t.insert(t.end(),atom_list.begin(),atom_list.end());
232 : } else {
233 0 : error(strings[i] + " is not a label plumed knows");
234 : }
235 : } else {
236 0 : error("atoms specified using @ symbol but no MOLINFO was available");
237 : }
238 : }
239 : }
240 : // here we check if the atom name is the name of a group
241 446265 : if(!ok) {
242 7623 : if(atoms.groups.count(strings[i])) {
243 : const auto m=atoms.groups.find(strings[i]);
244 404 : t.insert(t.end(),m->second.begin(),m->second.end());
245 : ok=true;
246 : }
247 : }
248 : // here we check if the atom name is the name of an added virtual atom
249 446624 : if(!ok) {
250 7219 : const ActionSet&actionSet(plumed.getActionSet());
251 6088236 : for(const auto & a : actionSet) {
252 6088236 : ActionWithVirtualAtom* c=dynamic_cast<ActionWithVirtualAtom*>(a.get());
253 6088236 : if(c)
254 6080098 : if(c->getLabel()==strings[i]) {
255 : ok=true;
256 7219 : t.push_back(c->getIndex());
257 : break;
258 : }
259 : }
260 : }
261 439809 : if(!ok) {
262 0 : error("it was not possible to interpret atom name " + strings[i]);
263 : }
264 : // plumed_massert(ok,"it was not possible to interpret atom name " + strings[i]);
265 : }
266 11914 : }
267 :
268 174709 : void ActionAtomistic::retrieveAtoms() {
269 174709 : pbc=atoms.pbc;
270 174709 : Colvar*cc=dynamic_cast<Colvar*>(this);
271 174709 : if(cc && cc->checkIsEnergy()) {
272 3989 : energy=atoms.getEnergy();
273 : }
274 174709 : if(donotretrieve) {
275 : return;
276 : }
277 168753 : chargesWereSet=atoms.chargesWereSet();
278 : const std::vector<Vector> & p(atoms.positions);
279 : const std::vector<double> & c(atoms.charges);
280 : const std::vector<double> & m(atoms.masses);
281 4759435 : for(unsigned j=0; j<indexes.size(); j++) {
282 4590682 : positions[j]=p[indexes[j].index()];
283 : }
284 4759435 : for(unsigned j=0; j<indexes.size(); j++) {
285 4590682 : charges[j]=c[indexes[j].index()];
286 : }
287 4759435 : for(unsigned j=0; j<indexes.size(); j++) {
288 4590682 : masses[j]=m[indexes[j].index()];
289 : }
290 : }
291 :
292 864 : void ActionAtomistic::setForcesOnAtoms(const std::vector<double>& forcesToApply, unsigned ind) {
293 864 : if(donotforce) {
294 : return;
295 : }
296 295850 : for(unsigned i=0; i<indexes.size(); ++i) {
297 294986 : forces[i][0]=forcesToApply[ind];
298 294986 : ind++;
299 294986 : forces[i][1]=forcesToApply[ind];
300 294986 : ind++;
301 294986 : forces[i][2]=forcesToApply[ind];
302 294986 : ind++;
303 : }
304 864 : virial(0,0)=forcesToApply[ind];
305 864 : ind++;
306 864 : virial(0,1)=forcesToApply[ind];
307 864 : ind++;
308 864 : virial(0,2)=forcesToApply[ind];
309 864 : ind++;
310 864 : virial(1,0)=forcesToApply[ind];
311 864 : ind++;
312 864 : virial(1,1)=forcesToApply[ind];
313 864 : ind++;
314 864 : virial(1,2)=forcesToApply[ind];
315 864 : ind++;
316 864 : virial(2,0)=forcesToApply[ind];
317 864 : ind++;
318 864 : virial(2,1)=forcesToApply[ind];
319 864 : ind++;
320 864 : virial(2,2)=forcesToApply[ind];
321 : plumed_dbg_assert( ind+1==forcesToApply.size());
322 : }
323 :
324 174246 : void ActionAtomistic::applyForces() {
325 174246 : if(donotforce) {
326 : return;
327 : }
328 168290 : std::vector<Vector>& f(atoms.forces);
329 168290 : Tensor& v(atoms.virial);
330 4739485 : for(unsigned j=0; j<indexes.size(); j++) {
331 4571195 : f[indexes[j].index()]+=forces[j];
332 : }
333 168290 : v+=virial;
334 168290 : atoms.forceOnEnergy+=forceOnEnergy;
335 168290 : if(extraCV.length()>0) {
336 48 : atoms.updateExtraCVForce(extraCV,forceOnExtraCV);
337 : }
338 : }
339 :
340 174709 : void ActionAtomistic::clearOutputForces() {
341 174709 : virial.zero();
342 174709 : if(donotforce) {
343 : return;
344 : }
345 168753 : Tools::set_to_zero(forces);
346 168753 : forceOnEnergy=0.0;
347 168753 : forceOnExtraCV=0.0;
348 : }
349 :
350 :
351 0 : void ActionAtomistic::readAtomsFromPDB(const PDB& pdb) {
352 0 : Colvar*cc=dynamic_cast<Colvar*>(this);
353 0 : if(cc && cc->checkIsEnergy()) {
354 0 : error("can't read energies from pdb files");
355 : }
356 :
357 0 : for(unsigned j=0; j<indexes.size(); j++) {
358 0 : if( indexes[j].index()>pdb.size() ) {
359 0 : error("there are not enough atoms in the input pdb file");
360 : }
361 0 : if( pdb.getAtomNumbers()[j].index()!=indexes[j].index() ) {
362 0 : error("there are atoms missing in the pdb file");
363 : }
364 0 : positions[j]=pdb.getPositions()[indexes[j].index()];
365 : }
366 0 : for(unsigned j=0; j<indexes.size(); j++) {
367 0 : charges[j]=pdb.getBeta()[indexes[j].index()];
368 : }
369 0 : for(unsigned j=0; j<indexes.size(); j++) {
370 0 : masses[j]=pdb.getOccupancy()[indexes[j].index()];
371 : }
372 0 : }
373 :
374 107745 : void ActionAtomistic::makeWhole() {
375 1160578 : for(unsigned j=0; j<positions.size()-1; ++j) {
376 : const Vector & first (positions[j]);
377 1052833 : Vector & second (positions[j+1]);
378 1052833 : second=first+pbcDistance(first,second);
379 : }
380 107745 : }
381 :
382 20022 : void ActionAtomistic::updateUniqueLocal() {
383 20022 : if(atoms.dd && atoms.shuffledAtoms>0) {
384 9706 : unique_local.clear();
385 91318 : for(auto pp=unique.begin(); pp!=unique.end(); ++pp) {
386 81612 : if(atoms.g2l[pp->index()]>=0) {
387 32944 : unique_local.push_back(*pp); // already sorted
388 : }
389 : }
390 : } else {
391 10316 : unique_local=unique; // copy
392 : }
393 20022 : }
394 :
395 : }
|
