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 10109 : ActionAtomistic::~ActionAtomistic() {
39 : // forget the pending request
40 10109 : atoms.remove(this);
41 10109 : }
42 :
43 10109 : ActionAtomistic::ActionAtomistic(const ActionOptions&ao):
44 : Action(ao),
45 10109 : lockRequestAtoms(false),
46 10109 : donotretrieve(false),
47 10109 : donotforce(false),
48 10109 : atoms(plumed.getAtoms())
49 : {
50 10109 : atoms.add(this);
51 : // if(atoms.getNatoms()==0) error("Cannot perform calculations involving atoms without atoms");
52 10109 : }
53 :
54 10199 : void ActionAtomistic::registerKeywords( Keywords& keys ) {
55 : (void) keys; // avoid warning
56 10199 : }
57 :
58 :
59 10145 : void ActionAtomistic::requestAtoms(const std::vector<AtomNumber> & a, const bool clearDep) {
60 10145 : plumed_massert(!lockRequestAtoms,"requested atom list can only be changed in the prepare() method");
61 10145 : int nat=a.size();
62 10145 : indexes=a;
63 10145 : positions.resize(nat);
64 10145 : forces.resize(nat);
65 10145 : masses.resize(nat);
66 10145 : charges.resize(nat);
67 10145 : int n=atoms.positions.size();
68 10145 : if(clearDep) clearDependencies();
69 : unique.clear();
70 1180256 : for(unsigned i=0; i<indexes.size(); i++) {
71 1170114 : if(indexes[i].index()>=n) { std::string num; Tools::convert( indexes[i].serial(),num ); error("atom " + num + " out of range"); }
72 1170111 : if(atoms.isVirtualAtom(indexes[i])) addDependency(atoms.getVirtualAtomsAction(indexes[i]));
73 : // only real atoms are requested to lower level Atoms class
74 1162890 : else unique.insert(indexes[i]);
75 : }
76 10144 : updateUniqueLocal();
77 10144 : atoms.unique.clear();
78 10144 : }
79 :
80 181662127 : Vector ActionAtomistic::pbcDistance(const Vector &v1,const Vector &v2)const {
81 181662127 : return pbc.distance(v1,v2);
82 : }
83 :
84 220567 : void ActionAtomistic::pbcApply(std::vector<Vector>& dlist, unsigned max_index)const {
85 220567 : pbc.apply(dlist, max_index);
86 220567 : }
87 :
88 195 : void ActionAtomistic::calculateNumericalDerivatives( ActionWithValue* a ) {
89 195 : calculateAtomicNumericalDerivatives( a, 0 );
90 195 : }
91 :
92 0 : void ActionAtomistic::changeBox( const Tensor& newbox ) {
93 0 : pbc.setBox( newbox );
94 0 : }
95 :
96 492 : void ActionAtomistic::calculateAtomicNumericalDerivatives( ActionWithValue* a, const unsigned& startnum ) {
97 492 : if(!a) {
98 279 : a=dynamic_cast<ActionWithValue*>(this);
99 279 : plumed_massert(a,"only Actions with a value can be differentiated");
100 : }
101 :
102 492 : const size_t nval=a->getNumberOfComponents();
103 492 : const size_t natoms=getNumberOfAtoms();
104 492 : std::vector<Vector> value(nval*natoms);
105 492 : std::vector<Tensor> valuebox(nval);
106 492 : std::vector<Vector> savedPositions(natoms);
107 : const double delta=std::sqrt(epsilon);
108 :
109 91856 : for(int i=0; i<natoms; i++) for(int k=0; k<3; k++) {
110 68523 : savedPositions[i][k]=positions[i][k];
111 68523 : positions[i][k]=positions[i][k]+delta;
112 68523 : a->calculate();
113 68523 : positions[i][k]=savedPositions[i][k];
114 202917 : for(int j=0; j<nval; j++) {
115 134394 : value[j*natoms+i][k]=a->getOutputQuantity(j);
116 : }
117 : }
118 492 : Tensor box(pbc.getBox());
119 6396 : for(int i=0; i<3; i++) for(int k=0; k<3; k++) {
120 4428 : double arg0=box(i,k);
121 209997 : for(int j=0; j<natoms; j++) positions[j]=pbc.realToScaled(positions[j]);
122 4428 : box(i,k)=box(i,k)+delta;
123 4428 : pbc.setBox(box);
124 209997 : for(int j=0; j<natoms; j++) positions[j]=pbc.scaledToReal(positions[j]);
125 4428 : a->calculate();
126 4428 : box(i,k)=arg0;
127 4428 : pbc.setBox(box);
128 209997 : for(int j=0; j<natoms; j++) positions[j]=savedPositions[j];
129 19512 : for(int j=0; j<nval; j++) valuebox[j](i,k)=a->getOutputQuantity(j);
130 : }
131 :
132 492 : a->calculate();
133 492 : a->clearDerivatives();
134 2168 : for(int j=0; j<nval; j++) {
135 1676 : Value* v=a->copyOutput(j);
136 : double ref=v->get();
137 1676 : if(v->hasDerivatives()) {
138 100117 : for(int i=0; i<natoms; i++) for(int k=0; k<3; k++) {
139 74532 : double d=(value[j*natoms+i][k]-ref)/delta;
140 74532 : v->addDerivative(startnum+3*i+k,d);
141 : }
142 741 : Tensor virial;
143 9633 : for(int i=0; i<3; i++) for(int k=0; k<3; k++)virial(i,k)= (valuebox[j](i,k)-ref)/delta;
144 : // BE CAREFUL WITH NON ORTHOROMBIC CELL
145 741 : virial=-matmul(box.transpose(),virial);
146 9633 : for(int i=0; i<3; i++) for(int k=0; k<3; k++) v->addDerivative(startnum+3*natoms+3*k+i,virial(k,i));
147 : }
148 : }
149 492 : }
150 :
151 11601 : void ActionAtomistic::parseAtomList(const std::string&key, std::vector<AtomNumber> &t) {
152 11601 : parseAtomList(key,-1,t);
153 11600 : }
154 :
155 13577 : void ActionAtomistic::parseAtomList(const std::string&key,const int num, std::vector<AtomNumber> &t) {
156 13577 : plumed_massert( keywords.style(key,"atoms") || keywords.style(key,"hidden"), "keyword " + key + " should be registered as atoms");
157 : std::vector<std::string> strings;
158 13577 : if( num<0 ) {
159 11601 : parseVector(key,strings);
160 11600 : if(strings.empty()) return;
161 : } else {
162 1976 : if ( !parseNumberedVector(key,num,strings) ) return;
163 : }
164 11075 : interpretAtomList( strings, t );
165 13577 : }
166 :
167 11274 : void ActionAtomistic::interpretAtomList(std::vector<std::string>& strings, std::vector<AtomNumber> &t) {
168 11274 : Tools::interpretRanges(strings); t.resize(0);
169 250699 : for(unsigned i=0; i<strings.size(); ++i) {
170 : AtomNumber atom;
171 239425 : bool ok=Tools::convertNoexcept(strings[i],atom); // this is converting strings to AtomNumbers
172 239425 : if(ok) t.push_back(atom);
173 : // here we check if this is a special symbol for MOLINFO
174 239425 : if( !ok && strings[i].compare(0,1,"@")==0 ) {
175 631 : std::string symbol=strings[i].substr(1);
176 631 : if(symbol=="allatoms") {
177 10 : const auto n=plumed.getAtoms().getNatoms() + plumed.getAtoms().getNVirtualAtoms();
178 10 : t.reserve(t.size()+n);
179 365 : for(unsigned i=0; i<n; i++) t.push_back(AtomNumber::index(i));
180 : ok=true;
181 621 : } else if(symbol=="mdatoms") {
182 4 : const auto n=plumed.getAtoms().getNatoms();
183 4 : t.reserve(t.size()+n);
184 228 : for(unsigned i=0; i<n; i++) t.push_back(AtomNumber::index(i));
185 : ok=true;
186 : } else {
187 617 : auto* moldat=plumed.getActionSet().selectLatest<GenericMolInfo*>(this);
188 617 : if( moldat ) {
189 617 : std::vector<AtomNumber> atom_list; moldat->interpretSymbol( symbol, atom_list );
190 617 : if( atom_list.size()>0 ) { ok=true; t.insert(t.end(),atom_list.begin(),atom_list.end()); }
191 0 : else { error(strings[i] + " is not a label plumed knows"); }
192 : } else {
193 0 : error("atoms specified using @ symbol but no MOLINFO was available");
194 : }
195 : }
196 : }
197 : // here we check if the atom name is the name of a group
198 238794 : if(!ok) {
199 7596 : if(atoms.groups.count(strings[i])) {
200 : const auto m=atoms.groups.find(strings[i]);
201 381 : t.insert(t.end(),m->second.begin(),m->second.end());
202 : ok=true;
203 : }
204 : }
205 : // here we check if the atom name is the name of an added virtual atom
206 239044 : if(!ok) {
207 7215 : const ActionSet&actionSet(plumed.getActionSet());
208 6088154 : for(const auto & a : actionSet) {
209 6088154 : ActionWithVirtualAtom* c=dynamic_cast<ActionWithVirtualAtom*>(a.get());
210 6088154 : if(c) if(c->getLabel()==strings[i]) {
211 : ok=true;
212 7215 : t.push_back(c->getIndex());
213 : break;
214 : }
215 : }
216 : }
217 232210 : if(!ok) error("it was not possible to interpret atom name " + strings[i]);
218 : // plumed_massert(ok,"it was not possible to interpret atom name " + strings[i]);
219 : }
220 11274 : }
221 :
222 158317 : void ActionAtomistic::retrieveAtoms() {
223 158317 : pbc=atoms.pbc;
224 158317 : Colvar*cc=dynamic_cast<Colvar*>(this);
225 158317 : if(cc && cc->checkIsEnergy()) energy=atoms.getEnergy();
226 158317 : if(donotretrieve) return;
227 154425 : chargesWereSet=atoms.chargesWereSet();
228 : const std::vector<Vector> & p(atoms.positions);
229 : const std::vector<double> & c(atoms.charges);
230 : const std::vector<double> & m(atoms.masses);
231 3688318 : for(unsigned j=0; j<indexes.size(); j++) positions[j]=p[indexes[j].index()];
232 3688318 : for(unsigned j=0; j<indexes.size(); j++) charges[j]=c[indexes[j].index()];
233 3688318 : for(unsigned j=0; j<indexes.size(); j++) masses[j]=m[indexes[j].index()];
234 : }
235 :
236 624 : void ActionAtomistic::setForcesOnAtoms(const std::vector<double>& forcesToApply, unsigned ind) {
237 624 : if(donotforce) return;
238 213170 : for(unsigned i=0; i<indexes.size(); ++i) {
239 212546 : forces[i][0]=forcesToApply[ind]; ind++;
240 212546 : forces[i][1]=forcesToApply[ind]; ind++;
241 212546 : forces[i][2]=forcesToApply[ind]; ind++;
242 : }
243 624 : virial(0,0)=forcesToApply[ind]; ind++;
244 624 : virial(0,1)=forcesToApply[ind]; ind++;
245 624 : virial(0,2)=forcesToApply[ind]; ind++;
246 624 : virial(1,0)=forcesToApply[ind]; ind++;
247 624 : virial(1,1)=forcesToApply[ind]; ind++;
248 624 : virial(1,2)=forcesToApply[ind]; ind++;
249 624 : virial(2,0)=forcesToApply[ind]; ind++;
250 624 : virial(2,1)=forcesToApply[ind]; ind++;
251 624 : virial(2,2)=forcesToApply[ind];
252 : plumed_dbg_assert( ind+1==forcesToApply.size());
253 : }
254 :
255 157858 : void ActionAtomistic::applyForces() {
256 157858 : if(donotforce) return;
257 153966 : std::vector<Vector>& f(atoms.forces);
258 153966 : Tensor& v(atoms.virial);
259 3668372 : for(unsigned j=0; j<indexes.size(); j++) f[indexes[j].index()]+=forces[j];
260 153966 : v+=virial;
261 153966 : atoms.forceOnEnergy+=forceOnEnergy;
262 153966 : if(extraCV.length()>0) atoms.updateExtraCVForce(extraCV,forceOnExtraCV);
263 : }
264 :
265 158317 : void ActionAtomistic::clearOutputForces() {
266 158317 : virial.zero();
267 158317 : if(donotforce) return;
268 3688318 : for(unsigned i=0; i<forces.size(); ++i) forces[i].zero();
269 154425 : forceOnEnergy=0.0;
270 154425 : forceOnExtraCV=0.0;
271 : }
272 :
273 :
274 0 : void ActionAtomistic::readAtomsFromPDB(const PDB& pdb) {
275 0 : Colvar*cc=dynamic_cast<Colvar*>(this);
276 0 : if(cc && cc->checkIsEnergy()) error("can't read energies from pdb files");
277 :
278 0 : for(unsigned j=0; j<indexes.size(); j++) {
279 0 : if( indexes[j].index()>pdb.size() ) error("there are not enough atoms in the input pdb file");
280 0 : if( pdb.getAtomNumbers()[j].index()!=indexes[j].index() ) error("there are atoms missing in the pdb file");
281 0 : positions[j]=pdb.getPositions()[indexes[j].index()];
282 : }
283 0 : for(unsigned j=0; j<indexes.size(); j++) charges[j]=pdb.getBeta()[indexes[j].index()];
284 0 : for(unsigned j=0; j<indexes.size(); j++) masses[j]=pdb.getOccupancy()[indexes[j].index()];
285 0 : }
286 :
287 106195 : void ActionAtomistic::makeWhole() {
288 1150846 : for(unsigned j=0; j<positions.size()-1; ++j) {
289 : const Vector & first (positions[j]);
290 1044651 : Vector & second (positions[j+1]);
291 1044651 : second=first+pbcDistance(first,second);
292 : }
293 106195 : }
294 :
295 19794 : void ActionAtomistic::updateUniqueLocal() {
296 : unique_local.clear();
297 19794 : if(atoms.dd && atoms.shuffledAtoms>0) {
298 91318 : for(auto pp=unique.begin(); pp!=unique.end(); ++pp) {
299 81612 : if(atoms.g2l[pp->index()]>=0) unique_local.insert(*pp);
300 : }
301 : } else {
302 10088 : unique_local.insert(unique.begin(),unique.end());
303 : }
304 19794 : }
305 :
306 : }
|
