Line data Source code
1 : /* +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
2 : Copyright (c) 2011-2019 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 "SetupMolInfo.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 : using namespace std;
37 :
38 : namespace PLMD {
39 :
40 2819 : ActionAtomistic::~ActionAtomistic() {
41 : // forget the pending request
42 2819 : atoms.remove(this);
43 2819 : }
44 :
45 2819 : ActionAtomistic::ActionAtomistic(const ActionOptions&ao):
46 : Action(ao),
47 : lockRequestAtoms(false),
48 : donotretrieve(false),
49 : donotforce(false),
50 14095 : atoms(plumed.getAtoms())
51 : {
52 2819 : atoms.add(this);
53 : // if(atoms.getNatoms()==0) error("Cannot perform calculations involving atoms without atoms");
54 2819 : }
55 :
56 2912 : void ActionAtomistic::registerKeywords( Keywords& keys ) {
57 : (void) keys; // avoid warning
58 2912 : }
59 :
60 :
61 2825 : void ActionAtomistic::requestAtoms(const vector<AtomNumber> & a, const bool clearDep) {
62 2825 : plumed_massert(!lockRequestAtoms,"requested atom list can only be changed in the prepare() method");
63 2825 : int nat=a.size();
64 2825 : indexes=a;
65 2825 : positions.resize(nat);
66 2825 : forces.resize(nat);
67 2825 : masses.resize(nat);
68 2825 : charges.resize(nat);
69 2825 : int n=atoms.positions.size();
70 2825 : if(clearDep) clearDependencies();
71 : unique.clear();
72 2224724 : for(unsigned i=0; i<indexes.size(); i++) {
73 1109542 : if(indexes[i].index()>=n) { std::string num; Tools::convert( indexes[i].serial(),num ); error("atom " + num + " out of range"); }
74 2219315 : if(atoms.isVirtualAtom(indexes[i])) addDependency(atoms.getVirtualAtomsAction(indexes[i]));
75 : // only real atoms are requested to lower level Atoms class
76 : else unique.insert(indexes[i]);
77 : }
78 2824 : updateUniqueLocal();
79 2824 : atoms.unique.clear();
80 2824 : }
81 :
82 182673070 : Vector ActionAtomistic::pbcDistance(const Vector &v1,const Vector &v2)const {
83 365338631 : return pbc.distance(v1,v2);
84 : }
85 :
86 240397 : void ActionAtomistic::pbcApply(std::vector<Vector>& dlist, unsigned max_index)const {
87 240397 : pbc.apply(dlist, max_index);
88 240413 : }
89 :
90 180 : void ActionAtomistic::calculateNumericalDerivatives( ActionWithValue* a ) {
91 180 : calculateAtomicNumericalDerivatives( a, 0 );
92 180 : }
93 :
94 0 : void ActionAtomistic::changeBox( const Tensor& newbox ) {
95 0 : pbc.setBox( newbox );
96 0 : }
97 :
98 472 : void ActionAtomistic::calculateAtomicNumericalDerivatives( ActionWithValue* a, const unsigned& startnum ) {
99 472 : if(!a) {
100 255 : a=dynamic_cast<ActionWithValue*>(this);
101 255 : plumed_massert(a,"only Actions with a value can be differentiated");
102 : }
103 :
104 : const int nval=a->getNumberOfComponents();
105 472 : const int natoms=getNumberOfAtoms();
106 472 : std::vector<Vector> value(nval*natoms);
107 472 : std::vector<Tensor> valuebox(nval);
108 472 : std::vector<Vector> savedPositions(natoms);
109 : const double delta=sqrt(epsilon);
110 :
111 63469 : for(int i=0; i<natoms; i++) for(int k=0; k<3; k++) {
112 188991 : savedPositions[i][k]=positions[i][k];
113 62997 : positions[i][k]=positions[i][k]+delta;
114 62997 : a->calculate();
115 125994 : positions[i][k]=savedPositions[i][k];
116 163503 : for(int j=0; j<nval; j++) {
117 201012 : value[j*natoms+i][k]=a->getOutputQuantity(j);
118 : }
119 : }
120 472 : Tensor box(pbc.getBox());
121 4720 : for(int i=0; i<3; i++) for(int k=0; k<3; k++) {
122 4248 : double arg0=box(i,k);
123 193239 : for(int j=0; j<natoms; j++) positions[j]=pbc.realToScaled(positions[j]);
124 4248 : box(i,k)=box(i,k)+delta;
125 4248 : pbc.setBox(box);
126 377982 : for(int j=0; j<natoms; j++) positions[j]=pbc.scaledToReal(positions[j]);
127 4248 : a->calculate();
128 4248 : box(i,k)=arg0;
129 4248 : pbc.setBox(box);
130 377982 : for(int j=0; j<natoms; j++) positions[j]=savedPositions[j];
131 28584 : for(int j=0; j<nval; j++) valuebox[j](i,k)=a->getOutputQuantity(j);
132 : }
133 :
134 472 : a->calculate();
135 472 : a->clearDerivatives();
136 1588 : for(int j=0; j<nval; j++) {
137 1588 : Value* v=a->copyOutput(j);
138 : double ref=v->get();
139 1588 : if(v->hasDerivatives()) {
140 69006 : for(int i=0; i<natoms; i++) for(int k=0; k<3; k++) {
141 138012 : double d=(value[j*natoms+i][k]-ref)/delta;
142 69006 : v->addDerivative(startnum+3*i+k,d);
143 : }
144 721 : Tensor virial;
145 7210 : for(int i=0; i<3; i++) for(int k=0; k<3; k++)virial(i,k)= (valuebox[j](i,k)-ref)/delta;
146 : // BE CAREFUL WITH NON ORTHOROMBIC CELL
147 721 : virial=-matmul(box.transpose(),virial);
148 7210 : 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));
149 : }
150 : }
151 472 : }
152 :
153 3843 : void ActionAtomistic::parseAtomList(const std::string&key, std::vector<AtomNumber> &t) {
154 3843 : parseAtomList(key,-1,t);
155 3842 : }
156 :
157 5781 : void ActionAtomistic::parseAtomList(const std::string&key,const int num, std::vector<AtomNumber> &t) {
158 11562 : plumed_massert( keywords.style(key,"atoms") || keywords.style(key,"hidden"), "keyword " + key + " should be registered as atoms");
159 : vector<string> strings;
160 5781 : if( num<0 ) {
161 3843 : parseVector(key,strings);
162 5837 : if(strings.empty()) return;
163 : } else {
164 1938 : if ( !parseNumberedVector(key,num,strings) ) return;
165 : }
166 3786 : interpretAtomList( strings, t );
167 : }
168 :
169 3974 : void ActionAtomistic::interpretAtomList( std::vector<std::string>& strings, std::vector<AtomNumber> &t) {
170 3974 : Tools::interpretRanges(strings); t.resize(0);
171 412662 : for(unsigned i=0; i<strings.size(); ++i) {
172 : AtomNumber atom;
173 202357 : bool ok=Tools::convert(strings[i],atom); // this is converting strings to AtomNumbers
174 202357 : if(ok) t.push_back(atom);
175 : // here we check if this is a special symbol for MOLINFO
176 203507 : if( !ok && strings[i].compare(0,1,"@")==0 ) {
177 567 : std::string symbol=strings[i].substr(1);
178 567 : if(symbol=="allatoms") {
179 4 : const auto n=plumed.getAtoms().getNatoms() + plumed.getAtoms().getNVirtualAtoms();
180 2 : t.reserve(t.size()+n);
181 440 : for(unsigned i=0; i<n; i++) t.push_back(AtomNumber::index(i));
182 : ok=true;
183 565 : } else if(symbol=="mdatoms") {
184 2 : const auto n=plumed.getAtoms().getNatoms();
185 2 : t.reserve(t.size()+n);
186 432 : for(unsigned i=0; i<n; i++) t.push_back(AtomNumber::index(i));
187 : ok=true;
188 : } else {
189 1126 : vector<SetupMolInfo*> moldat=plumed.getActionSet().select<SetupMolInfo*>();
190 563 : if( moldat.size()>0 ) {
191 563 : vector<AtomNumber> atom_list; moldat[0]->interpretSymbol( symbol, atom_list );
192 563 : ok=true; t.insert(t.end(),atom_list.begin(),atom_list.end());
193 : } else {
194 0 : error("atoms specified using @ symbol but no MOLINFO was available");
195 : }
196 : }
197 : }
198 : // here we check if the atom name is the name of a group
199 202357 : if(!ok) {
200 1166 : if(atoms.groups.count(strings[i])) {
201 348 : const auto m=atoms.groups.find(strings[i]);
202 348 : t.insert(t.end(),m->second.begin(),m->second.end());
203 : ok=true;
204 : }
205 : }
206 : // here we check if the atom name is the name of an added virtual atom
207 202357 : if(!ok) {
208 235 : const ActionSet&actionSet(plumed.getActionSet());
209 1813 : for(const auto & a : actionSet) {
210 1578 : ActionWithVirtualAtom* c=dynamic_cast<ActionWithVirtualAtom*>(a.get());
211 2068 : if(c) if(c->getLabel()==strings[i]) {
212 : ok=true;
213 470 : t.push_back(c->getIndex());
214 235 : break;
215 : }
216 : }
217 : }
218 202357 : if(!ok) error("it was not possible to interpret atom name " + strings[i]);
219 : // plumed_massert(ok,"it was not possible to interpret atom name " + strings[i]);
220 : }
221 3974 : }
222 :
223 :
224 140149 : void ActionAtomistic::retrieveAtoms() {
225 140149 : pbc=atoms.pbc;
226 140149 : Colvar*cc=dynamic_cast<Colvar*>(this);
227 245710 : if(cc && cc->checkIsEnergy()) energy=atoms.getEnergy();
228 280298 : if(donotretrieve) return;
229 273162 : chargesWereSet=atoms.chargesWereSet();
230 : const vector<Vector> & p(atoms.positions);
231 : const vector<double> & c(atoms.charges);
232 : const vector<double> & m(atoms.masses);
233 9137868 : for(unsigned j=0; j<indexes.size(); j++) positions[j]=p[indexes[j].index()];
234 9001287 : for(unsigned j=0; j<indexes.size(); j++) charges[j]=c[indexes[j].index()];
235 9001287 : for(unsigned j=0; j<indexes.size(); j++) masses[j]=m[indexes[j].index()];
236 : }
237 :
238 624 : void ActionAtomistic::setForcesOnAtoms( const std::vector<double>& forcesToApply, unsigned ind ) {
239 1248 : if(donotforce) return;
240 425716 : for(unsigned i=0; i<indexes.size(); ++i) {
241 425092 : forces[i][0]=forcesToApply[ind]; ind++;
242 425092 : forces[i][1]=forcesToApply[ind]; ind++;
243 425092 : forces[i][2]=forcesToApply[ind]; ind++;
244 : }
245 1248 : virial(0,0)=forcesToApply[ind]; ind++;
246 1248 : virial(0,1)=forcesToApply[ind]; ind++;
247 1248 : virial(0,2)=forcesToApply[ind]; ind++;
248 1248 : virial(1,0)=forcesToApply[ind]; ind++;
249 1248 : virial(1,1)=forcesToApply[ind]; ind++;
250 1248 : virial(1,2)=forcesToApply[ind]; ind++;
251 1248 : virial(2,0)=forcesToApply[ind]; ind++;
252 1248 : virial(2,1)=forcesToApply[ind]; ind++;
253 1248 : virial(2,2)=forcesToApply[ind];
254 : plumed_dbg_assert( ind+1==forcesToApply.size());
255 : }
256 :
257 139720 : void ActionAtomistic::applyForces() {
258 279440 : if(donotforce) return;
259 136152 : vector<Vector> & f(atoms.forces);
260 136152 : Tensor & v(atoms.virial);
261 9078609 : for(unsigned j=0; j<indexes.size(); j++) f[indexes[j].index()]+=forces[j];
262 136152 : v+=virial;
263 136152 : atoms.forceOnEnergy+=forceOnEnergy;
264 136152 : if(extraCV.length()>0) atoms.updateExtraCVForce(extraCV,forceOnExtraCV);
265 : }
266 :
267 140149 : void ActionAtomistic::clearOutputForces() {
268 140149 : virial.zero();
269 280298 : if(donotforce) return;
270 6046385 : for(unsigned i=0; i<forces.size(); ++i)forces[i].zero();
271 136581 : forceOnEnergy=0.0;
272 136581 : forceOnExtraCV=0.0;
273 : }
274 :
275 :
276 0 : void ActionAtomistic::readAtomsFromPDB( const PDB& pdb ) {
277 0 : Colvar*cc=dynamic_cast<Colvar*>(this);
278 0 : if(cc && cc->checkIsEnergy()) error("can't read energies from pdb files");
279 :
280 0 : for(unsigned j=0; j<indexes.size(); j++) {
281 0 : if( indexes[j].index()>pdb.size() ) error("there are not enough atoms in the input pdb file");
282 0 : if( pdb.getAtomNumbers()[j].index()!=indexes[j].index() ) error("there are atoms missing in the pdb file");
283 0 : positions[j]=pdb.getPositions()[indexes[j].index()];
284 : }
285 0 : for(unsigned j=0; j<indexes.size(); j++) charges[j]=pdb.getBeta()[indexes[j].index()];
286 0 : for(unsigned j=0; j<indexes.size(); j++) masses[j]=pdb.getOccupancy()[indexes[j].index()];
287 0 : }
288 :
289 93508 : void ActionAtomistic::makeWhole() {
290 1269922 : for(unsigned j=0; j<positions.size()-1; ++j) {
291 : const Vector & first (positions[j]);
292 494699 : Vector & second (positions[j+1]);
293 494699 : second=first+pbcDistance(first,second);
294 : }
295 93508 : }
296 :
297 12380 : void ActionAtomistic::updateUniqueLocal() {
298 : unique_local.clear();
299 24760 : if(atoms.dd && atoms.shuffledAtoms>0) {
300 94640 : for(auto pp=unique.begin(); pp!=unique.end(); ++pp) {
301 150960 : if(atoms.g2l[pp->index()]>=0) unique_local.insert(*pp);
302 : }
303 : } else {
304 : unique_local.insert(unique.begin(),unique.end());
305 : }
306 12380 : }
307 :
308 5874 : }
|
