Line data Source code
1 : /* +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
2 : Copyright (c) 2014-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 "core/ActionAtomistic.h"
23 : #include "core/ActionPilot.h"
24 : #include "core/ActionRegister.h"
25 : #include "tools/Vector.h"
26 : #include "tools/AtomNumber.h"
27 : #include "tools/Tools.h"
28 : #include "tools/Pbc.h"
29 : #include "core/PlumedMain.h"
30 : #include "core/ActionSet.h"
31 : #include "core/GenericMolInfo.h"
32 :
33 : #include <vector>
34 :
35 : namespace PLMD {
36 : namespace generic {
37 :
38 : //+PLUMEDOC GENERIC WRAPAROUND
39 : /*
40 : Rebuild periodic boundary conditions around chosen atoms.
41 :
42 : This action modifies the position of the atoms indicated by ATOMS by shifting them by lattice vectors so that they are
43 : as close as possible to the atoms indicated by AROUND. More precisely, for every atom i
44 : in the ATOMS list the following procedure is performed:
45 : - The atom j among those in the AROUND list is searched that is closest to atom i.
46 : - The atom i is replaced with its periodic image that is closest to atom j.
47 :
48 : This action works similarly to [WHOLEMOLECULES](WHOLEMOLECULES.md) in that it replaces atoms coordinate. Notice that only
49 : atoms specified with ATOMS are replaced, and that, at variance with [WHOLEMOLECULES](WHOLEMOLECULES.md),
50 : the order in which atoms are specified is irrelevant.
51 :
52 : This is often convenient at a post processing stage (using the driver), but sometime
53 : it is required during the simulation if collective variables need atoms to be in a specific periodic image.
54 :
55 : !!! caution "modifies stored positions
56 :
57 : This directive modifies the stored position at the precise moment it is executed. This means that only collective variables which are below it in
58 : the input script will see the corrected positions. As a general rule, put it at the top of the input file. Also, unless you know exactly what you are doing,
59 : leave the default stride (1), so that this action is performed at every MD step.
60 :
61 : The computational cost of this action grows with the product
62 : of the size of the two lists (ATOMS and AROUND), so this action can become very expensive.
63 : If you are using it to analyze a trajectory this is usually not a big problem. If you use it to
64 : analyze a simulation on the fly, e.g. with [DUMPATOMS](DUMPATOMS.md) to store a properly wrapped trajectory,
65 : consider using the STRIDE keyword here (with great care).
66 :
67 :
68 : ## Examples
69 :
70 : This command instructs plumed to move all the ions to their periodic image that is as close as possible to
71 : the rna group.
72 :
73 : ```plumed
74 : rna: GROUP ATOMS=1-100
75 : ions: GROUP ATOMS=101-110
76 : # first make the rna molecule whole
77 : WHOLEMOLECULES ENTITY0=rna
78 : WRAPAROUND ATOMS=ions AROUND=rna
79 : DUMPATOMS FILE=dump.xyz ATOMS=rna,ions
80 : ```
81 :
82 : In case you want to do it during a simulation and you only care about wrapping the ions in
83 : the `dump.xyz` file, you can use the following input:
84 :
85 : ```plumed
86 : # add some restraint that do not require molecules to be whole:
87 : a: TORSION ATOMS=1,2,10,11
88 : RESTRAINT ARG=a AT=0.0 KAPPA=5
89 :
90 :
91 : # then do the things that are required for dumping the trajectory
92 : # notice that they are all done every 100 steps, so as not to
93 : # unnecessarily overload the calculation
94 :
95 : rna: GROUP ATOMS=1-100
96 : ions: GROUP ATOMS=101-110
97 : # first make the rna molecule whole
98 : WHOLEMOLECULES ENTITY0=rna STRIDE=100
99 : WRAPAROUND ATOMS=ions AROUND=rna STRIDE=100
100 : DUMPATOMS FILE=dump.xyz ATOMS=rna,ions STRIDE=100
101 : ```
102 :
103 : Notice that if the biased variable requires a molecule to be whole, you might have to put
104 : the [WHOLEMOLECULES](WHOLEMOLECULES.md) command before computing that variable and leave the default STRIDE=1.
105 :
106 : This command instructs plumed to center all atoms around the center of mass of a solute molecule.
107 :
108 : ```plumed
109 : solute: GROUP ATOMS=1-100
110 : all: GROUP ATOMS=1-1000
111 : # center of the solute:
112 : # notice that since plumed 2.2 this also works if the
113 : # solute molecule is broken
114 : com: COM ATOMS=solute
115 : # notice that we wrap around a single atom. this should be fast
116 : WRAPAROUND ATOMS=all AROUND=com
117 : DUMPATOMS FILE=dump.xyz ATOMS=all
118 : ```
119 :
120 : Notice that whereas [WHOLEMOLECULES](WHOLEMOLECULES.md) is designed to make molecules whole,
121 : WRAPAROUND can easily break molecules. In the last example,
122 : if solvent (atoms 101-1000) is made e.g. of water, then water
123 : molecules could be broken by WRAPAROUND (hydrogen could end up
124 : in an image and oxygen in another one).
125 : One solution is to use [WHOLEMOLECULES](WHOLEMOLECULES.md) on _all_ the water molecules
126 : after WRAPAROUND. This is tedious. A better solution is to use the
127 : GROUPBY option which is going
128 : to consider the atoms listed in ATOMS as a list of groups
129 : each of size GROUPBY. The first atom of the group will be brought
130 : close to the AROUND atoms. The following atoms of the group
131 : will be just brought close to the first atom of the group.
132 : Assuming that oxygen is the first atom of each water molecules,
133 : in the following examples all the water oxygen atoms will be brought
134 : close to the solute, and all the hydrogen atoms will be kept close
135 : to their related oxygen.
136 :
137 : ```plumed
138 : solute: GROUP ATOMS=1-100
139 : water: GROUP ATOMS=101-1000
140 : com: COM ATOMS=solute
141 : # notice that we wrap around a single atom. this should be fast
142 : WRAPAROUND ATOMS=solute AROUND=com
143 : # notice that we wrap around a single atom. this should be fast
144 : WRAPAROUND ATOMS=water AROUND=com GROUPBY=3
145 : DUMPATOMS FILE=dump.xyz ATOMS=solute,water
146 : ```
147 :
148 : */
149 : //+ENDPLUMEDOC
150 :
151 :
152 : class WrapAround:
153 : public ActionPilot,
154 : public ActionAtomistic {
155 : // cppcheck-suppress duplInheritedMember
156 : std::vector<Vector> refatoms;
157 : // distances from the atom being wrapped to all the reference atoms
158 : std::vector<Vector> distances;
159 : std::vector<std::pair<std::size_t,std::size_t> > p_atoms;
160 : std::vector<std::pair<std::size_t,std::size_t> > p_reference;
161 : unsigned groupby;
162 : bool pair_;
163 : public:
164 : explicit WrapAround(const ActionOptions&ao);
165 : static void registerKeywords( Keywords& keys );
166 584 : bool actionHasForces() override {
167 584 : return false;
168 : }
169 : void calculate() override;
170 590 : void apply() override {}
171 : };
172 :
173 : PLUMED_REGISTER_ACTION(WrapAround,"WRAPAROUND")
174 :
175 11 : void WrapAround::registerKeywords( Keywords& keys ) {
176 11 : Action::registerKeywords( keys );
177 11 : ActionAtomistic::registerKeywords( keys );
178 11 : ActionPilot::registerKeywords( keys );
179 11 : keys.add("compulsory","STRIDE","1","the frequency with which molecules are reassembled. Unless you are completely certain about what you are doing leave this set equal to 1!");
180 11 : keys.add("atoms","AROUND","reference atoms");
181 11 : keys.add("atoms","ATOMS","wrapped atoms");
182 11 : keys.add("compulsory","GROUPBY","1","group atoms so as not to break molecules");
183 11 : keys.addFlag("PAIR", false, "Pair atoms in AROUND and ATOMS groups");
184 11 : }
185 :
186 7 : WrapAround::WrapAround(const ActionOptions&ao):
187 : Action(ao),
188 : ActionPilot(ao),
189 : ActionAtomistic(ao),
190 7 : groupby(1),
191 7 : pair_(false) {
192 7 : std::vector<AtomNumber> atoms;
193 7 : parseAtomList("ATOMS",atoms);
194 7 : std::vector<AtomNumber> reference;
195 14 : parseAtomList("AROUND",reference);
196 7 : if(reference.empty()) {
197 0 : error("no reference atoms given with AROUND");
198 : }
199 7 : parse("GROUPBY",groupby);
200 7 : if(groupby<1) {
201 1 : error("GROUPBY should be at least 1");
202 : }
203 6 : parseFlag("PAIR", pair_);
204 :
205 6 : log.printf(" atoms in reference :");
206 13 : for(unsigned j=0; j<reference.size(); ++j) {
207 7 : log.printf(" %d",reference[j].serial() );
208 : }
209 6 : log.printf("\n");
210 6 : log.printf(" atoms to be wrapped :");
211 422 : for(unsigned j=0; j<atoms.size(); ++j) {
212 416 : log.printf(" %d",atoms[j].serial() );
213 : }
214 6 : log.printf("\n");
215 6 : if(groupby>1) {
216 1 : log<<" atoms will be grouped by "<<groupby<<"\n";
217 : }
218 6 : if(pair_) {
219 0 : log.printf(" pairing atoms and references\n");
220 : }
221 :
222 6 : if(atoms.size()%groupby!=0) {
223 0 : error("number of atoms should be a multiple of groupby option");
224 : }
225 : // additional checks with PAIR
226 6 : if(pair_ && atoms.size()!=reference.size()*groupby) {
227 0 : error("with PAIR you must have: #ATOMS = #AROUND * #GROUPBY");
228 : }
229 :
230 6 : checkRead();
231 :
232 : // do not remove duplicates with pair
233 6 : if(!pair_) {
234 6 : if(groupby<=1) {
235 5 : Tools::removeDuplicates(atoms);
236 : }
237 6 : Tools::removeDuplicates(reference);
238 : }
239 :
240 6 : std::vector<AtomNumber> merged(atoms.size()+reference.size());
241 6 : merge(atoms.begin(),atoms.end(),reference.begin(),reference.end(),merged.begin());
242 6 : p_atoms.resize( atoms.size() );
243 422 : for(unsigned i=0; i<atoms.size(); ++i) {
244 416 : p_atoms[i] = getValueIndices( atoms[i] );
245 : }
246 6 : refatoms.resize( reference.size() );
247 6 : p_reference.resize( reference.size() );
248 13 : for(unsigned i=0; i<reference.size(); ++i) {
249 7 : p_reference[i] = getValueIndices( reference[i] );
250 : }
251 6 : Tools::removeDuplicates(merged);
252 6 : requestAtoms(merged);
253 : doNotRetrieve();
254 : doNotForce();
255 14 : }
256 :
257 590 : void WrapAround::calculate() {
258 590 : const unsigned nref=p_reference.size();
259 1185 : for(unsigned j=0; j<nref; ++j) {
260 595 : refatoms[j] = getGlobalPosition(p_reference[j]);
261 : }
262 590 : if(!pair_) {
263 590 : distances.resize(nref);
264 : }
265 : const Pbc& pbc=getPbc();
266 : #ifdef __PLUMED_PBC_WHILE
267 : const bool orthorombic=false;
268 : #else
269 590 : const bool orthorombic=pbc.isOrthorombic();
270 : #endif
271 590 : const Vector diag(pbc.getBox()(0,0),pbc.getBox()(1,1),pbc.getBox()(2,2));
272 590 : const Vector invdiag(pbc.getInvBox()(0,0),pbc.getInvBox()(1,1),pbc.getInvBox()(2,2));
273 15265 : for(unsigned i=0; i<p_atoms.size(); i+=groupby) {
274 14675 : Vector second, first=getGlobalPosition(p_atoms[i]);
275 : int closest=-1;
276 14675 : if(pair_) {
277 0 : closest = i/groupby;
278 : } else {
279 : double mindist2=std::numeric_limits<double>::max();
280 14675 : if(orthorombic) {
281 : // same operations as Pbc::distance, written out so that the loop is inlined
282 29900 : for(unsigned j=0; j<nref; ++j) {
283 : double distance2=0.0;
284 60900 : for(unsigned k=0; k<3; ++k) {
285 45675 : const double d=Tools::pbc((refatoms[j][k]-first[k])*invdiag[k])*diag[k];
286 45675 : distance2+=d*d;
287 : }
288 15225 : if(distance2<mindist2) {
289 : mindist2=distance2;
290 14955 : closest=j;
291 : }
292 : }
293 : } else {
294 : // all the distances to the reference atoms are computed in one call,
295 : // the result is the same as calling pbcDistance(first,refatoms[j]) for each j
296 0 : for(unsigned j=0; j<nref; ++j) {
297 0 : distances[j]=delta(first,refatoms[j]);
298 : }
299 0 : pbc.apply(VectorView(&distances[0][0],nref));
300 0 : for(unsigned j=0; j<nref; ++j) {
301 0 : const double distance2=modulo2(distances[j]);
302 0 : if(distance2<mindist2) {
303 : mindist2=distance2;
304 0 : closest=j;
305 : }
306 : }
307 : }
308 14675 : plumed_massert(closest>=0,"closest not found");
309 : }
310 14675 : second=refatoms[closest];
311 : // place first atom of the group
312 14675 : first=second+pbcDistance(second,first);
313 14675 : setGlobalPosition(p_atoms[i],first);
314 : // then place other atoms close to the first of the group
315 15170 : for(unsigned j=1; j<groupby; j++) {
316 495 : second=getGlobalPosition(p_atoms[i+j]);
317 495 : setGlobalPosition( p_atoms[i+j], first+pbcDistance(first,second) );
318 : }
319 : }
320 590 : }
321 :
322 :
323 :
324 : }
325 :
326 : }
|
