1 /* +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
2 Copyright (c) 2014-2021 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 "core/Atoms.h"
29 #include "core/PlumedMain.h"
30 #include "core/ActionSet.h"
31 #include "core/GenericMolInfo.h"
32
33 #include <vector>
34 #include <string>
35 #include <limits>
36
37 using namespace std;
38
39 namespace PLMD {
40 namespace generic {
41
42 //+PLUMEDOC GENERIC WRAPAROUND
43 /*
44 Rebuild periodic boundary conditions around chosen atoms.
45
46
47 Modify position of atoms indicated by ATOMS by shifting them by lattice vectors so that they are
48 as close as possible to the atoms indicated by AROUND. More precisely, for every atom i
49 in the ATOMS list the following procedure is performed:
50 - The atom j among those in the AROUND list is searched that is closest to atom i.
51 - The atom i is replaced with its periodic image that is closest to atom j.
52
53 This action works similarly to \ref WHOLEMOLECULES in that it replaces atoms coordinate. Notice that only
54 atoms specified with ATOMS are replaced, and that, at variance with \ref WHOLEMOLECULES,
55 the order in which atoms are specified is irrelevant.
56
57 This is often convenient at a post processing stage (using the \ref driver), but sometime
58 it is required during the simulation if collective variables need atoms to be in a specific periodic image.
59
60 \attention This directive modifies the stored position at the precise moment it is executed. This means that only collective variables which are below it in 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, leave the default stride (1), so that this action is performed at every MD step.
61
62 Consider that the computational cost grows with the product
63 of the size of the two lists (ATOMS and AROUND), so that this action can become very expensive.
64 If you are using it to analyze a trajectory this is usually not a big problem. If you use it to
65 analyze a simulation on the fly, e.g. with \ref DUMPATOMS to store a properly wrapped trajectory,
66 consider the possibility of using the STRIDE keyword here (with great care).
67 \par Examples
68
69 This command instructs plumed to move all the ions to their periodic image that is as close as possible to
70 the rna group.
71
72 \plumedfile
73 rna: GROUP ATOMS=1-100
74 ions: GROUP ATOMS=101-110
75 # first make the rna molecule whole
76 WHOLEMOLECULES ENTITY0=rna
77 WRAPAROUND ATOMS=ions AROUND=rna
78 DUMPATOMS FILE=dump.xyz ATOMS=rna,ions
79 \endplumedfile
80
81 In case you want to do it during a simulation and you only care about wrapping the ions in
82 the `dump.xyz` file, you can use the following:
83
84 \plumedfile
85 # add some restraint that do not require molecules to be whole:
86 a: TORSION ATOMS=1,2,10,11
87 RESTRAINT ARG=a AT=0.0 KAPPA=5
88
89
90 # then do the things that are required for dumping the trajectory
91 # notice that they are all done every 100 steps, so as not to
92 # unnecessarily overload the calculation
93
94 rna: GROUP ATOMS=1-100
95 ions: GROUP ATOMS=101-110
96 # first make the rna molecule whole
97 WHOLEMOLECULES ENTITY0=rna STRIDE=100
98 WRAPAROUND ATOMS=ions AROUND=rna STRIDE=100
99 DUMPATOMS FILE=dump.xyz ATOMS=rna,ions STRIDE=100
100 \endplumedfile
101
102 Notice that if the biased variable requires a molecule to be whole, you might have to put
103 just the \ref WHOLEMOLECULES command before computing that variable and leave the default STRIDE=1.
104
105 This command instructs plumed to center all atoms around the center of mass of a solute molecule.
106
107 \plumedfile
108 solute: GROUP ATOMS=1-100
109 all: GROUP ATOMS=1-1000
110 # center of the solute:
111 # notice that since plumed 2.2 this also works if the
112 # solute molecule is broken
113 com: COM ATOMS=solute
114 # notice that we wrap around a single atom. this should be fast
115 WRAPAROUND ATOMS=all AROUND=com
116 DUMPATOMS FILE=dump.xyz ATOMS=all
117 \endplumedfile
118
119 Notice that whereas \ref WHOLEMOLECULES is designed to make molecules whole,
120 \ref WRAPAROUND can easily break molecules. In the last example,
121 if solvent (atoms 101-1000) is made e.g. of water, then water
122 molecules could be broken by \ref WRAPAROUND (hydrogen could end up
123 in an image and oxygen in another one).
124 One solution is to use \ref WHOLEMOLECULES on _all_ the water molecules
125 after \ref WRAPAROUND. This is tedious. A better solution is to use the
126 GROUPBY option which is going
127 to consider the atoms listed in ATOMS as a list of groups
128 each of size GROUPBY. The first atom of the group will be brought
129 close to the AROUND atoms. The following atoms of the group
130 will be just brought close to the first atom of the group.
131 Assuming that oxygen is the first atom of each water molecules,
132 in the following examples all the water oxygen atoms will be brought
133 close to the solute, and all the hydrogen atoms will be kept close
134 to their related oxygen.
135
136 \plumedfile
137 solute: GROUP ATOMS=1-100
138 water: GROUP ATOMS=101-1000
139 com: COM ATOMS=solute
140 # notice that we wrap around a single atom. this should be fast
141 WRAPAROUND ATOMS=solute AROUND=com
142 # notice that we wrap around a single atom. this should be fast
143 WRAPAROUND ATOMS=water AROUND=com GROUPBY=3
144 DUMPATOMS FILE=dump.xyz ATOMS=solute,water
145 \endplumedfile
146
147 */
148 //+ENDPLUMEDOC
149
150
151 class WrapAround:
152 public ActionPilot,
153 public ActionAtomistic
154 {
155 vector<AtomNumber> atoms;
156 vector<AtomNumber> reference;
157 unsigned groupby;
158 public:
159 explicit WrapAround(const ActionOptions&ao);
160 static void registerKeywords( Keywords& keys );
161 void calculate() override;
apply()162 void apply() override {}
163 };
164
165 PLUMED_REGISTER_ACTION(WrapAround,"WRAPAROUND")
166
registerKeywords(Keywords & keys)167 void WrapAround::registerKeywords( Keywords& keys ) {
168 Action::registerKeywords( keys );
169 ActionAtomistic::registerKeywords( keys );
170 ActionPilot::registerKeywords( keys );
171 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!");
172 keys.add("atoms","AROUND","reference atoms");
173 keys.add("atoms","ATOMS","wrapped atoms");
174 keys.add("compulsory","GROUPBY","1","group atoms so as not to break molecules");
175 }
176
WrapAround(const ActionOptions & ao)177 WrapAround::WrapAround(const ActionOptions&ao):
178 Action(ao),
179 ActionPilot(ao),
180 ActionAtomistic(ao),
181 groupby(1)
182 {
183 parseAtomList("ATOMS",atoms);
184 parseAtomList("AROUND",reference);
185 parse("GROUPBY",groupby);
186
187 log.printf(" atoms in reference :");
188 for(unsigned j=0; j<reference.size(); ++j) log.printf(" %d",reference[j].serial() );
189 log.printf("\n");
190 log.printf(" atoms to be wrapped :");
191 for(unsigned j=0; j<atoms.size(); ++j) log.printf(" %d",atoms[j].serial() );
192 log.printf("\n");
193 if(groupby>1) log<<" atoms will be grouped by "<<groupby<<"\n";
194
195 if(atoms.size()%groupby!=0) error("number of atoms should be a multiple of groupby option");
196
197 checkRead();
198
199 if(groupby<=1) Tools::removeDuplicates(atoms);
200 Tools::removeDuplicates(reference);
201
202 vector<AtomNumber> merged(atoms.size()+reference.size());
203 merge(atoms.begin(),atoms.end(),reference.begin(),reference.end(),merged.begin());
204 Tools::removeDuplicates(merged);
205 requestAtoms(merged);
206 doNotRetrieve();
207 doNotForce();
208 }
209
calculate()210 void WrapAround::calculate() {
211 for(unsigned i=0; i<atoms.size(); i+=groupby) {
212 Vector & first (modifyGlobalPosition(atoms[i]));
213 double mindist2=std::numeric_limits<double>::max();
214 int closest=-1;
215 for(unsigned j=0; j<reference.size(); ++j) {
216 Vector & second (modifyGlobalPosition(reference[j]));
217 Vector distance=pbcDistance(first,second);
218 double distance2=modulo2(distance);
219 if(distance2<mindist2) {
220 mindist2=distance2;
221 closest=j;
222 }
223 }
224 plumed_massert(closest>=0,"closest not found");
225 Vector & second (modifyGlobalPosition(reference[closest]));
226 // place first atom of the group
227 first=second+pbcDistance(second,first);
228 // then place other atoms close to the first of the group
229 for(unsigned j=1; j<groupby; j++) {
230 Vector & second (modifyGlobalPosition(atoms[i+j]));
231 second=first+pbcDistance(first,second);
232 }
233 }
234 }
235
236
237
238 }
239
240 }
241