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 2310 : ActionAtomistic::~ActionAtomistic() {
41 : // forget the pending request
42 2310 : atoms.remove(this);
43 2310 : delete(&pbc);
44 2310 : }
45 :
46 2310 : ActionAtomistic::ActionAtomistic(const ActionOptions&ao):
47 : Action(ao),
48 2310 : pbc(*new(Pbc)),
49 : lockRequestAtoms(false),
50 : donotretrieve(false),
51 : donotforce(false),
52 9240 : atoms(plumed.getAtoms())
53 : {
54 2310 : atoms.add(this);
55 : // if(atoms.getNatoms()==0) error("Cannot perform calculations involving atoms without atoms");
56 2310 : }
57 :
58 2393 : void ActionAtomistic::registerKeywords( Keywords& keys ) {
59 : (void) keys; // avoid warning
60 2393 : }
61 :
62 :
63 2238 : void ActionAtomistic::requestAtoms(const vector<AtomNumber> & a) {
64 2238 : plumed_massert(!lockRequestAtoms,"requested atom list can only be changed in the prepare() method");
65 2238 : int nat=a.size();
66 2238 : indexes=a;
67 2238 : positions.resize(nat);
68 2238 : forces.resize(nat);
69 2238 : masses.resize(nat);
70 2238 : charges.resize(nat);
71 2238 : int n=atoms.positions.size();
72 2238 : clearDependencies();
73 : unique.clear();
74 1835916 : for(unsigned i=0; i<indexes.size(); i++) {
75 610482 : if(indexes[i].index()>=n) error("atom out of range");
76 1221144 : if(atoms.isVirtualAtom(indexes[i])) addDependency(atoms.getVirtualAtomsAction(indexes[i]));
77 : // only real atoms are requested to lower level Atoms class
78 : else unique.insert(indexes[i]);
79 : }
80 2237 : updateUniqueLocal();
81 2237 : atoms.unique.clear();
82 2237 : }
83 :
84 127141753 : Vector ActionAtomistic::pbcDistance(const Vector &v1,const Vector &v2)const {
85 254283506 : return pbc.distance(v1,v2);
86 : }
87 :
88 234812 : void ActionAtomistic::pbcApply(std::vector<Vector>& dlist, unsigned max_index)const {
89 234812 : pbc.apply(dlist, max_index);
90 234812 : }
91 :
92 175 : void ActionAtomistic::calculateNumericalDerivatives( ActionWithValue* a ) {
93 175 : calculateAtomicNumericalDerivatives( a, 0 );
94 175 : }
95 :
96 0 : void ActionAtomistic::changeBox( const Tensor& newbox ) {
97 0 : pbc.setBox( newbox );
98 0 : }
99 :
100 467 : void ActionAtomistic::calculateAtomicNumericalDerivatives( ActionWithValue* a, const unsigned& startnum ) {
101 467 : if(!a) {
102 250 : a=dynamic_cast<ActionWithValue*>(this);
103 250 : plumed_massert(a,"only Actions with a value can be differentiated");
104 : }
105 :
106 : const int nval=a->getNumberOfComponents();
107 467 : const int natoms=getNumberOfAtoms();
108 467 : std::vector<Vector> value(nval*natoms);
109 467 : std::vector<Tensor> valuebox(nval);
110 467 : std::vector<Vector> savedPositions(natoms);
111 : const double delta=sqrt(epsilon);
112 :
113 66255 : for(int i=0; i<natoms; i++) for(int k=0; k<3; k++) {
114 148023 : savedPositions[i][k]=positions[i][k];
115 49341 : positions[i][k]=positions[i][k]+delta;
116 49341 : a->calculate();
117 98682 : positions[i][k]=savedPositions[i][k];
118 204981 : for(int j=0; j<nval; j++) {
119 155640 : value[j*natoms+i][k]=a->getOutputQuantity(j);
120 : }
121 : }
122 467 : Tensor box(pbc.getBox());
123 6071 : for(int i=0; i<3; i++) for(int k=0; k<3; k++) {
124 4203 : double arg0=box(i,k);
125 300249 : for(int j=0; j<natoms; j++) positions[j]=pbc.realToScaled(positions[j]);
126 4203 : box(i,k)=box(i,k)+delta;
127 4203 : pbc.setBox(box);
128 300249 : for(int j=0; j<natoms; j++) positions[j]=pbc.scaledToReal(positions[j]);
129 4203 : a->calculate();
130 4203 : box(i,k)=arg0;
131 4203 : pbc.setBox(box);
132 300249 : for(int j=0; j<natoms; j++) positions[j]=savedPositions[j];
133 32607 : for(int j=0; j<nval; j++) valuebox[j](i,k)=a->getOutputQuantity(j);
134 : }
135 :
136 467 : a->calculate();
137 467 : a->clearDerivatives();
138 3623 : for(int j=0; j<nval; j++) {
139 1578 : Value* v=a->copyOutput(j);
140 : double ref=v->get();
141 1578 : if(v->hasDerivatives()) {
142 76925 : for(int i=0; i<natoms; i++) for(int k=0; k<3; k++) {
143 114312 : double d=(value[j*natoms+i][k]-ref)/delta;
144 57156 : v->addDerivative(startnum+3*i+k,d);
145 : }
146 717 : Tensor virial;
147 9321 : for(int i=0; i<3; i++) for(int k=0; k<3; k++)virial(i,k)= (valuebox[j](i,k)-ref)/delta;
148 : // BE CAREFUL WITH NON ORTHOROMBIC CELL
149 717 : virial=-matmul(box.transpose(),virial);
150 9321 : 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));
151 : }
152 : }
153 467 : }
154 :
155 3222 : void ActionAtomistic::parseAtomList(const std::string&key, std::vector<AtomNumber> &t) {
156 3222 : parseAtomList(key,-1,t);
157 3221 : }
158 :
159 4817 : void ActionAtomistic::parseAtomList(const std::string&key,const int num, std::vector<AtomNumber> &t) {
160 14451 : plumed_massert( keywords.style(key,"atoms") || keywords.style(key,"hidden"), "keyword " + key + " should be registered as atoms");
161 3146 : vector<string> strings;
162 4817 : if( num<0 ) {
163 3222 : parseVector(key,strings);
164 4892 : if(strings.empty()) return;
165 : } else {
166 1595 : if ( !parseNumberedVector(key,num,strings) ) return;
167 : }
168 3145 : interpretAtomList( strings, t );
169 : }
170 :
171 3323 : void ActionAtomistic::interpretAtomList( std::vector<std::string>& strings, std::vector<AtomNumber> &t) {
172 3323 : Tools::interpretRanges(strings); t.resize(0);
173 479221 : for(unsigned i=0; i<strings.size(); ++i) {
174 : AtomNumber atom;
175 157525 : bool ok=Tools::convert(strings[i],atom); // this is converting strings to AtomNumbers
176 157525 : if(ok) t.push_back(atom);
177 : // here we check if this is a special symbol for MOLINFO
178 158507 : if( !ok && strings[i].compare(0,1,"@")==0 ) {
179 502 : std::size_t dot=strings[i].find_first_of("@"); std::string symbol=strings[i].substr(dot+1);
180 1004 : vector<SetupMolInfo*> moldat=plumed.getActionSet().select<SetupMolInfo*>();
181 502 : if( moldat.size()>0 ) {
182 502 : vector<AtomNumber> atom_list; moldat[0]->interpretSymbol( symbol, atom_list );
183 1004 : if( atom_list.size()>0 ) { ok=true; t.insert(t.end(),atom_list.begin(),atom_list.end()); }
184 0 : else { error(strings[i] + " is not a label plumed knows"); }
185 : } else {
186 0 : error("atoms specified using @ symbol but no MOLINFO was available");
187 : }
188 : }
189 : // here we check if the atom name is the name of a group
190 157525 : if(!ok) {
191 480 : if(atoms.groups.count(strings[i])) {
192 298 : const auto m=atoms.groups.find(strings[i]);
193 298 : t.insert(t.end(),m->second.begin(),m->second.end());
194 : ok=true;
195 : }
196 : }
197 : // here we check if the atom name is the name of an added virtual atom
198 157227 : if(!ok) {
199 182 : const ActionSet&actionSet(plumed.getActionSet());
200 1175 : for(const auto & a : actionSet) {
201 1175 : ActionWithVirtualAtom* c=dynamic_cast<ActionWithVirtualAtom*>(a);
202 1795 : if(c) if(c->getLabel()==strings[i]) {
203 : ok=true;
204 364 : t.push_back(c->getIndex());
205 : break;
206 : }
207 : }
208 : }
209 157343 : if(!ok) error("it was not possible to interpret atom name " + strings[i]);
210 : // plumed_massert(ok,"it was not possible to interpret atom name " + strings[i]);
211 : }
212 3323 : }
213 :
214 :
215 92770 : void ActionAtomistic::retrieveAtoms() {
216 92770 : pbc=atoms.pbc;
217 92770 : Colvar*cc=dynamic_cast<Colvar*>(this);
218 170401 : if(cc && cc->checkIsEnergy()) energy=atoms.getEnergy();
219 92770 : if(donotretrieve) return;
220 183370 : chargesWereSet=atoms.chargesWereSet();
221 : const vector<Vector> & p(atoms.positions);
222 : const vector<double> & c(atoms.charges);
223 : const vector<double> & m(atoms.masses);
224 7666814 : for(unsigned j=0; j<indexes.size(); j++) positions[j]=p[indexes[j].index()];
225 7666814 : for(unsigned j=0; j<indexes.size(); j++) charges[j]=c[indexes[j].index()];
226 7666814 : for(unsigned j=0; j<indexes.size(); j++) masses[j]=m[indexes[j].index()];
227 : }
228 :
229 599 : void ActionAtomistic::setForcesOnAtoms( const std::vector<double>& forcesToApply, unsigned ind ) {
230 599 : if(donotforce) return;
231 638389 : for(unsigned i=0; i<indexes.size(); ++i) {
232 424794 : forces[i][0]=forcesToApply[ind]; ind++;
233 424794 : forces[i][1]=forcesToApply[ind]; ind++;
234 424794 : forces[i][2]=forcesToApply[ind]; ind++;
235 : }
236 1198 : virial(0,0)=forcesToApply[ind]; ind++;
237 1198 : virial(0,1)=forcesToApply[ind]; ind++;
238 1198 : virial(0,2)=forcesToApply[ind]; ind++;
239 1198 : virial(1,0)=forcesToApply[ind]; ind++;
240 1198 : virial(1,1)=forcesToApply[ind]; ind++;
241 1198 : virial(1,2)=forcesToApply[ind]; ind++;
242 1198 : virial(2,0)=forcesToApply[ind]; ind++;
243 1198 : virial(2,1)=forcesToApply[ind]; ind++;
244 1198 : virial(2,2)=forcesToApply[ind];
245 : plumed_dbg_assert( ind+1==forcesToApply.size());
246 : }
247 :
248 92371 : void ActionAtomistic::applyForces() {
249 92371 : if(donotforce) return;
250 91382 : vector<Vector> & f(atoms.forces);
251 91382 : Tensor & v(atoms.virial);
252 7598444 : for(unsigned j=0; j<indexes.size(); j++) f[indexes[j].index()]+=forces[j];
253 91382 : v+=virial;
254 91382 : atoms.forceOnEnergy+=forceOnEnergy;
255 : }
256 :
257 92770 : void ActionAtomistic::clearOutputForces() {
258 92770 : virial.zero();
259 92770 : if(donotforce) return;
260 5797153 : for(unsigned i=0; i<forces.size(); ++i)forces[i].zero();
261 91781 : forceOnEnergy=0.0;
262 : }
263 :
264 :
265 0 : void ActionAtomistic::readAtomsFromPDB( const PDB& pdb ) {
266 0 : Colvar*cc=dynamic_cast<Colvar*>(this);
267 0 : if(cc && cc->checkIsEnergy()) error("can't read energies from pdb files");
268 :
269 0 : for(unsigned j=0; j<indexes.size(); j++) {
270 0 : if( indexes[j].index()>pdb.size() ) error("there are not enough atoms in the input pdb file");
271 0 : if( pdb.getAtomNumbers()[j].index()!=indexes[j].index() ) error("there are atoms missing in the pdb file");
272 0 : positions[j]=pdb.getPositions()[indexes[j].index()];
273 : }
274 0 : for(unsigned j=0; j<indexes.size(); j++) charges[j]=pdb.getBeta()[indexes[j].index()];
275 0 : for(unsigned j=0; j<indexes.size(); j++) masses[j]=pdb.getOccupancy()[indexes[j].index()];
276 0 : }
277 :
278 52109 : void ActionAtomistic::makeWhole() {
279 1017967 : for(unsigned j=0; j<positions.size()-1; ++j) {
280 : const Vector & first (positions[j]);
281 304583 : Vector & second (positions[j+1]);
282 304583 : second=first+pbcDistance(first,second);
283 : }
284 52109 : }
285 :
286 11625 : void ActionAtomistic::updateUniqueLocal() {
287 : unique_local.clear();
288 23250 : if(atoms.dd && atoms.shuffledAtoms>0) {
289 79984 : for(auto pp=unique.begin(); pp!=unique.end(); ++pp) {
290 140952 : if(atoms.g2l[pp->index()]>=0) unique_local.insert(*pp);
291 : }
292 : } else {
293 : unique_local.insert(unique.begin(),unique.end());
294 : }
295 11625 : }
296 :
297 4839 : }
|
