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 "core/PlumedMain.h"
29 : #include "core/ActionSet.h"
30 : #include "core/GenericMolInfo.h"
31 :
32 : #include <vector>
33 :
34 : namespace PLMD {
35 : namespace generic {
36 :
37 : //+PLUMEDOC GENERIC WRAPAROUND
38 : /*
39 : Rebuild periodic boundary conditions around chosen atoms.
40 :
41 : This action modifies the position of the atoms indicated by ATOMS by shifting them by lattice vectors so that they are
42 : as close as possible to the atoms indicated by AROUND. More precisely, for every atom i
43 : in the ATOMS list the following procedure is performed:
44 : - The atom j among those in the AROUND list is searched that is closest to atom i.
45 : - The atom i is replaced with its periodic image that is closest to atom j.
46 :
47 : This action works similarly to [WHOLEMOLECULES](WHOLEMOLECULES.md) in that it replaces atoms coordinate. Notice that only
48 : atoms specified with ATOMS are replaced, and that, at variance with [WHOLEMOLECULES](WHOLEMOLECULES.md),
49 : the order in which atoms are specified is irrelevant.
50 :
51 : This is often convenient at a post processing stage (using the driver), but sometime
52 : it is required during the simulation if collective variables need atoms to be in a specific periodic image.
53 :
54 : !!! caution "modifies stored positions
55 :
56 : This directive modifies the stored position at the precise moment it is executed. This means that only collective variables which are below it in
57 : 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,
58 : leave the default stride (1), so that this action is performed at every MD step.
59 :
60 : The computational cost of this action grows with the product
61 : of the size of the two lists (ATOMS and AROUND), so this action can become very expensive.
62 : If you are using it to analyze a trajectory this is usually not a big problem. If you use it to
63 : analyze a simulation on the fly, e.g. with [DUMPATOMS](DUMPATOMS.md) to store a properly wrapped trajectory,
64 : consider using the STRIDE keyword here (with great care).
65 :
66 :
67 : ## 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 : ```plumed
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 : ```
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 input:
83 :
84 : ```plumed
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 : ```
101 :
102 : Notice that if the biased variable requires a molecule to be whole, you might have to put
103 : the [WHOLEMOLECULES](WHOLEMOLECULES.md) 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 : ```plumed
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 : ```
118 :
119 : Notice that whereas [WHOLEMOLECULES](WHOLEMOLECULES.md) is designed to make molecules whole,
120 : 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 WRAPAROUND (hydrogen could end up
123 : in an image and oxygen in another one).
124 : One solution is to use [WHOLEMOLECULES](WHOLEMOLECULES.md) on _all_ the water molecules
125 : after 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 : ```plumed
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 : ```
146 :
147 : */
148 : //+ENDPLUMEDOC
149 :
150 :
151 : class WrapAround:
152 : public ActionPilot,
153 : public ActionAtomistic {
154 : // cppcheck-suppress duplInheritedMember
155 : std::vector<Vector> refatoms;
156 : std::vector<std::pair<std::size_t,std::size_t> > p_atoms;
157 : std::vector<std::pair<std::size_t,std::size_t> > p_reference;
158 : unsigned groupby;
159 : bool pair_;
160 : public:
161 : explicit WrapAround(const ActionOptions&ao);
162 : static void registerKeywords( Keywords& keys );
163 0 : bool actionHasForces() override {
164 0 : return false;
165 : }
166 : void calculate() override;
167 590 : void apply() override {}
168 : };
169 :
170 : PLUMED_REGISTER_ACTION(WrapAround,"WRAPAROUND")
171 :
172 8 : void WrapAround::registerKeywords( Keywords& keys ) {
173 8 : Action::registerKeywords( keys );
174 8 : ActionAtomistic::registerKeywords( keys );
175 8 : ActionPilot::registerKeywords( keys );
176 8 : 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!");
177 8 : keys.add("atoms","AROUND","reference atoms");
178 8 : keys.add("atoms","ATOMS","wrapped atoms");
179 8 : keys.add("compulsory","GROUPBY","1","group atoms so as not to break molecules");
180 8 : keys.addFlag("PAIR", false, "Pair atoms in AROUND and ATOMS groups");
181 8 : }
182 :
183 6 : WrapAround::WrapAround(const ActionOptions&ao):
184 : Action(ao),
185 : ActionPilot(ao),
186 : ActionAtomistic(ao),
187 6 : groupby(1),
188 6 : pair_(false) {
189 : std::vector<AtomNumber> atoms;
190 12 : parseAtomList("ATOMS",atoms);
191 : std::vector<AtomNumber> reference;
192 6 : parseAtomList("AROUND",reference);
193 6 : parse("GROUPBY",groupby);
194 6 : parseFlag("PAIR", pair_);
195 :
196 6 : log.printf(" atoms in reference :");
197 13 : for(unsigned j=0; j<reference.size(); ++j) {
198 7 : log.printf(" %d",reference[j].serial() );
199 : }
200 6 : log.printf("\n");
201 6 : log.printf(" atoms to be wrapped :");
202 422 : for(unsigned j=0; j<atoms.size(); ++j) {
203 416 : log.printf(" %d",atoms[j].serial() );
204 : }
205 6 : log.printf("\n");
206 6 : if(groupby>1) {
207 1 : log<<" atoms will be grouped by "<<groupby<<"\n";
208 : }
209 6 : if(pair_) {
210 0 : log.printf(" pairing atoms and references\n");
211 : }
212 :
213 6 : if(atoms.size()%groupby!=0) {
214 0 : error("number of atoms should be a multiple of groupby option");
215 : }
216 : // additional checks with PAIR
217 6 : if(pair_ && atoms.size()!=reference.size()*groupby) {
218 0 : error("with PAIR you must have: #ATOMS = #AROUND * #GROUPBY");
219 : }
220 :
221 6 : checkRead();
222 :
223 : // do not remove duplicates with pair
224 6 : if(!pair_) {
225 6 : if(groupby<=1) {
226 5 : Tools::removeDuplicates(atoms);
227 : }
228 6 : Tools::removeDuplicates(reference);
229 : }
230 :
231 6 : std::vector<AtomNumber> merged(atoms.size()+reference.size());
232 6 : merge(atoms.begin(),atoms.end(),reference.begin(),reference.end(),merged.begin());
233 6 : p_atoms.resize( atoms.size() );
234 422 : for(unsigned i=0; i<atoms.size(); ++i) {
235 416 : p_atoms[i] = getValueIndices( atoms[i] );
236 : }
237 6 : refatoms.resize( reference.size() );
238 6 : p_reference.resize( reference.size() );
239 13 : for(unsigned i=0; i<reference.size(); ++i) {
240 7 : p_reference[i] = getValueIndices( reference[i] );
241 : }
242 6 : Tools::removeDuplicates(merged);
243 6 : requestAtoms(merged);
244 : doNotRetrieve();
245 : doNotForce();
246 6 : }
247 :
248 590 : void WrapAround::calculate() {
249 1185 : for(unsigned j=0; j<p_reference.size(); ++j) {
250 595 : refatoms[j] = getGlobalPosition(p_reference[j]);
251 : }
252 :
253 15265 : for(unsigned i=0; i<p_atoms.size(); i+=groupby) {
254 14675 : Vector second, first=getGlobalPosition(p_atoms[i]);
255 : double mindist2=std::numeric_limits<double>::max();
256 : int closest=-1;
257 14675 : if(pair_) {
258 0 : closest = i/groupby;
259 : } else {
260 29900 : for(unsigned j=0; j<p_reference.size(); ++j) {
261 15225 : second=refatoms[j];
262 : const Vector distance=pbcDistance(first,second);
263 : const double distance2=modulo2(distance);
264 15225 : if(distance2<mindist2) {
265 : mindist2=distance2;
266 14955 : closest=j;
267 : }
268 : }
269 14675 : plumed_massert(closest>=0,"closest not found");
270 : }
271 14675 : second=refatoms[closest];
272 : // place first atom of the group
273 14675 : first=second+pbcDistance(second,first);
274 14675 : setGlobalPosition(p_atoms[i],first);
275 : // then place other atoms close to the first of the group
276 15170 : for(unsigned j=1; j<groupby; j++) {
277 495 : second=getGlobalPosition(p_atoms[i+j]);
278 495 : setGlobalPosition( p_atoms[i+j], first+pbcDistance(first,second) );
279 : }
280 : }
281 590 : }
282 :
283 :
284 :
285 : }
286 :
287 : }
|
