Line data Source code
1 : /* +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
2 : Copyright (c) 2011-2018 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 2712 : ActionAtomistic::~ActionAtomistic() {
41 : // forget the pending request
42 1356 : atoms.remove(this);
43 1356 : delete(&pbc);
44 1356 : }
45 :
46 1356 : ActionAtomistic::ActionAtomistic(const ActionOptions&ao):
47 : Action(ao),
48 1356 : pbc(*new(Pbc)),
49 : lockRequestAtoms(false),
50 : donotretrieve(false),
51 : donotforce(false),
52 2712 : atoms(plumed.getAtoms())
53 : {
54 1356 : atoms.add(this);
55 : // if(atoms.getNatoms()==0) error("Cannot perform calculations involving atoms without atoms");
56 1356 : }
57 :
58 1450 : void ActionAtomistic::registerKeywords( Keywords& keys ) {
59 : (void) keys; // avoid warning
60 1450 : }
61 :
62 :
63 1303 : void ActionAtomistic::requestAtoms(const vector<AtomNumber> & a) {
64 1303 : plumed_massert(!lockRequestAtoms,"requested atom list can only be changed in the prepare() method");
65 1303 : int nat=a.size();
66 1303 : indexes=a;
67 1303 : positions.resize(nat);
68 1303 : forces.resize(nat);
69 1303 : masses.resize(nat);
70 1303 : charges.resize(nat);
71 1303 : int n=atoms.positions.size();
72 1303 : clearDependencies();
73 1303 : unique.clear();
74 153846 : for(unsigned i=0; i<indexes.size(); i++) {
75 152543 : if(indexes[i].index()>=n) error("atom out of range");
76 152543 : if(atoms.isVirtualAtom(indexes[i])) addDependency(atoms.getVirtualAtomsAction(indexes[i]));
77 : // only real atoms are requested to lower level Atoms class
78 152409 : else unique.insert(indexes[i]);
79 : }
80 :
81 1303 : }
82 :
83 8581984 : Vector ActionAtomistic::pbcDistance(const Vector &v1,const Vector &v2)const {
84 8581984 : return pbc.distance(v1,v2);
85 : }
86 :
87 191411 : void ActionAtomistic::pbcApply(std::vector<Vector>& dlist, unsigned max_index)const {
88 191411 : pbc.apply(dlist, max_index);
89 191436 : }
90 :
91 186 : void ActionAtomistic::calculateNumericalDerivatives( ActionWithValue* a ) {
92 186 : calculateAtomicNumericalDerivatives( a, 0 );
93 186 : }
94 :
95 0 : void ActionAtomistic::changeBox( const Tensor& newbox ) {
96 0 : pbc.setBox( newbox );
97 0 : }
98 :
99 403 : void ActionAtomistic::calculateAtomicNumericalDerivatives( ActionWithValue* a, const unsigned& startnum ) {
100 403 : if(!a) {
101 186 : a=dynamic_cast<ActionWithValue*>(this);
102 186 : plumed_massert(a,"only Actions with a value can be differentiated");
103 : }
104 :
105 403 : const int nval=a->getNumberOfComponents();
106 403 : const int natoms=getNumberOfAtoms();
107 403 : std::vector<Vector> value(nval*natoms);
108 806 : std::vector<Tensor> valuebox(nval);
109 806 : std::vector<Vector> savedPositions(natoms);
110 403 : const double delta=sqrt(epsilon);
111 :
112 46900 : for(int i=0; i<natoms; i++) for(int k=0; k<3; k++) {
113 46497 : savedPositions[i][k]=positions[i][k];
114 46497 : positions[i][k]=positions[i][k]+delta;
115 46497 : a->calculate();
116 46497 : positions[i][k]=savedPositions[i][k];
117 101091 : for(int j=0; j<nval; j++) {
118 54594 : value[j*natoms+i][k]=a->getOutputQuantity(j);
119 : }
120 : }
121 403 : Tensor box(pbc.getBox());
122 4030 : for(int i=0; i<3; i++) for(int k=0; k<3; k++) {
123 3627 : double arg0=box(i,k);
124 3627 : for(int j=0; j<natoms; j++) positions[j]=pbc.realToScaled(positions[j]);
125 3627 : box(i,k)=box(i,k)+delta;
126 3627 : pbc.setBox(box);
127 3627 : for(int j=0; j<natoms; j++) positions[j]=pbc.scaledToReal(positions[j]);
128 3627 : a->calculate();
129 3627 : box(i,k)=arg0;
130 3627 : pbc.setBox(box);
131 3627 : for(int j=0; j<natoms; j++) positions[j]=savedPositions[j];
132 3627 : for(int j=0; j<nval; j++) valuebox[j](i,k)=a->getOutputQuantity(j);
133 : }
134 :
135 403 : a->calculate();
136 403 : a->clearDerivatives();
137 923 : for(int j=0; j<nval; j++) {
138 520 : Value* v=a->copyOutput(j);
139 520 : double ref=v->get();
140 520 : if(v->hasDerivatives()) {
141 55114 : for(int i=0; i<natoms; i++) for(int k=0; k<3; k++) {
142 54594 : double d=(value[j*natoms+i][k]-ref)/delta;
143 54594 : v->addDerivative(startnum+3*i+k,d);
144 : }
145 520 : Tensor virial;
146 520 : for(int i=0; i<3; i++) for(int k=0; k<3; k++)virial(i,k)= (valuebox[j](i,k)-ref)/delta;
147 : // BE CAREFUL WITH NON ORTHOROMBIC CELL
148 520 : virial=-matmul(box.transpose(),virial);
149 520 : 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));
150 : }
151 403 : }
152 403 : }
153 :
154 1276 : void ActionAtomistic::parseAtomList(const std::string&key, std::vector<AtomNumber> &t) {
155 1276 : parseAtomList(key,-1,t);
156 1276 : }
157 :
158 2704 : void ActionAtomistic::parseAtomList(const std::string&key,const int num, std::vector<AtomNumber> &t) {
159 2704 : plumed_massert( keywords.style(key,"atoms") || keywords.style(key,"hidden"), "keyword " + key + " should be registered as atoms");
160 2704 : vector<string> strings;
161 2704 : if( num<0 ) {
162 1276 : parseVector(key,strings);
163 1276 : if(strings.empty()) return;
164 : } else {
165 1428 : if ( !parseNumberedVector(key,num,strings) ) return;
166 : }
167 2271 : interpretAtomList( strings, t );
168 : }
169 :
170 2410 : void ActionAtomistic::interpretAtomList( std::vector<std::string>& strings, std::vector<AtomNumber> &t) {
171 2410 : Tools::interpretRanges(strings); t.resize(0);
172 74386 : for(unsigned i=0; i<strings.size(); ++i) {
173 71976 : AtomNumber atom;
174 71976 : bool ok=Tools::convert(strings[i],atom); // this is converting strings to AtomNumbers
175 71976 : if(ok) t.push_back(atom);
176 : // here we check if this is a special symbol for MOLINFO
177 71976 : if( !ok && strings[i].compare(0,1,"@")==0 ) {
178 364 : std::size_t dot=strings[i].find_first_of("@"); std::string symbol=strings[i].substr(dot+1);
179 728 : vector<SetupMolInfo*> moldat=plumed.getActionSet().select<SetupMolInfo*>();
180 364 : if( moldat.size()>0 ) {
181 364 : vector<AtomNumber> atom_list; moldat[0]->interpretSymbol( symbol, atom_list );
182 364 : if( atom_list.size()>0 ) { ok=true; t.insert(t.end(),atom_list.begin(),atom_list.end()); }
183 0 : else { error(strings[i] + " is not a label plumed knows"); }
184 : } else {
185 0 : error("atoms specified using @ symbol but no MOLINFO was available");
186 364 : }
187 : }
188 : // here we check if the atom name is the name of a group
189 71976 : if(!ok) {
190 368 : if(atoms.groups.count(strings[i])) {
191 236 : map<string,vector<AtomNumber> >::const_iterator m=atoms.groups.find(strings[i]);
192 236 : t.insert(t.end(),m->second.begin(),m->second.end());
193 236 : ok=true;
194 : }
195 : }
196 : // here we check if the atom name is the name of an added virtual atom
197 71976 : if(!ok) {
198 132 : const ActionSet&actionSet(plumed.getActionSet());
199 573 : for(ActionSet::const_iterator a=actionSet.begin(); a!=actionSet.end(); ++a) {
200 573 : ActionWithVirtualAtom* c=dynamic_cast<ActionWithVirtualAtom*>(*a);
201 573 : if(c) if(c->getLabel()==strings[i]) {
202 132 : ok=true;
203 132 : t.push_back(c->getIndex());
204 132 : break;
205 : }
206 : }
207 : }
208 71976 : if(!ok) error("it was not possible to interpret atom name " + strings[i]);
209 : // plumed_massert(ok,"it was not possible to interpret atom name " + strings[i]);
210 : }
211 2410 : }
212 :
213 :
214 55928 : void ActionAtomistic::retrieveAtoms() {
215 55928 : pbc=atoms.pbc;
216 55928 : Colvar*cc=dynamic_cast<Colvar*>(this);
217 55928 : if(cc && cc->checkIsEnergy()) energy=atoms.getEnergy();
218 111856 : if(donotretrieve) return;
219 55658 : chargesWereSet=atoms.chargesWereSet();
220 55658 : const vector<Vector> & p(atoms.positions);
221 55658 : const vector<double> & c(atoms.charges);
222 55658 : const vector<double> & m(atoms.masses);
223 55658 : for(unsigned j=0; j<indexes.size(); j++) positions[j]=p[indexes[j].index()];
224 55658 : for(unsigned j=0; j<indexes.size(); j++) charges[j]=c[indexes[j].index()];
225 55658 : for(unsigned j=0; j<indexes.size(); j++) masses[j]=m[indexes[j].index()];
226 : }
227 :
228 140 : void ActionAtomistic::setForcesOnAtoms( const std::vector<double>& forcesToApply, unsigned ind ) {
229 280 : if(donotforce) return;
230 78541 : for(unsigned i=0; i<indexes.size(); ++i) {
231 78401 : forces[i][0]=forcesToApply[ind]; ind++;
232 78401 : forces[i][1]=forcesToApply[ind]; ind++;
233 78401 : forces[i][2]=forcesToApply[ind]; ind++;
234 : }
235 140 : virial(0,0)=forcesToApply[ind]; ind++;
236 140 : virial(0,1)=forcesToApply[ind]; ind++;
237 140 : virial(0,2)=forcesToApply[ind]; ind++;
238 140 : virial(1,0)=forcesToApply[ind]; ind++;
239 140 : virial(1,1)=forcesToApply[ind]; ind++;
240 140 : virial(1,2)=forcesToApply[ind]; ind++;
241 140 : virial(2,0)=forcesToApply[ind]; ind++;
242 140 : virial(2,1)=forcesToApply[ind]; ind++;
243 140 : virial(2,2)=forcesToApply[ind];
244 : plumed_dbg_assert( ind+1==forcesToApply.size());
245 : }
246 :
247 55604 : void ActionAtomistic::applyForces() {
248 111208 : if(donotforce) return;
249 55430 : vector<Vector> & f(atoms.forces);
250 55430 : Tensor & v(atoms.virial);
251 55430 : for(unsigned j=0; j<indexes.size(); j++) f[indexes[j].index()]+=forces[j];
252 55430 : v+=virial;
253 55430 : atoms.forceOnEnergy+=forceOnEnergy;
254 : }
255 :
256 55928 : void ActionAtomistic::clearOutputForces() {
257 55928 : virial.zero();
258 111856 : if(donotforce) return;
259 55754 : for(unsigned i=0; i<forces.size(); ++i)forces[i].zero();
260 55754 : forceOnEnergy=0.0;
261 : }
262 :
263 :
264 0 : void ActionAtomistic::readAtomsFromPDB( const PDB& pdb ) {
265 0 : Colvar*cc=dynamic_cast<Colvar*>(this);
266 0 : if(cc && cc->checkIsEnergy()) error("can't read energies from pdb files");
267 :
268 0 : for(unsigned j=0; j<indexes.size(); j++) {
269 0 : if( indexes[j].index()>pdb.size() ) error("there are not enough atoms in the input pdb file");
270 0 : if( pdb.getAtomNumbers()[j].index()!=indexes[j].index() ) error("there are atoms missing in the pdb file");
271 0 : positions[j]=pdb.getPositions()[indexes[j].index()];
272 : }
273 0 : for(unsigned j=0; j<indexes.size(); j++) charges[j]=pdb.getBeta()[indexes[j].index()];
274 0 : for(unsigned j=0; j<indexes.size(); j++) masses[j]=pdb.getOccupancy()[indexes[j].index()];
275 0 : }
276 :
277 32112 : void ActionAtomistic::makeWhole() {
278 194898 : for(unsigned j=0; j<positions.size()-1; ++j) {
279 162786 : const Vector & first (positions[j]);
280 162786 : Vector & second (positions[j+1]);
281 162786 : second=first+pbcDistance(first,second);
282 : }
283 32112 : }
284 :
285 2523 : }
|
