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]
55 : > 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.
56 :
57 : The computational cost of this action grows with the product
58 : of the size of the two lists (ATOMS and AROUND), so this action can become very expensive.
59 : If you are using it to analyze a trajectory this is usually not a big problem. If you use it to
60 : analyze a simulation on the fly, e.g. with [DUMPATOMS](DUMPATOMS.md) to store a properly wrapped trajectory,
61 : consider using the STRIDE keyword here (with great care).
62 :
63 :
64 : ## Examples
65 :
66 : This command instructs plumed to move all the ions to their periodic image that is as close as possible to
67 : the rna group.
68 :
69 : ```plumed
70 : rna: GROUP ATOMS=1-100
71 : ions: GROUP ATOMS=101-110
72 : # first make the rna molecule whole
73 : WHOLEMOLECULES ENTITY0=rna
74 : WRAPAROUND ATOMS=ions AROUND=rna
75 : DUMPATOMS FILE=dump.xyz ATOMS=rna,ions
76 : ```
77 :
78 : In case you want to do it during a simulation and you only care about wrapping the ions in
79 : the `dump.xyz` file, you can use the following input:
80 :
81 : ```plumed
82 : # add some restraint that do not require molecules to be whole:
83 : a: TORSION ATOMS=1,2,10,11
84 : RESTRAINT ARG=a AT=0.0 KAPPA=5
85 :
86 :
87 : # then do the things that are required for dumping the trajectory
88 : # notice that they are all done every 100 steps, so as not to
89 : # unnecessarily overload the calculation
90 :
91 : rna: GROUP ATOMS=1-100
92 : ions: GROUP ATOMS=101-110
93 : # first make the rna molecule whole
94 : WHOLEMOLECULES ENTITY0=rna STRIDE=100
95 : WRAPAROUND ATOMS=ions AROUND=rna STRIDE=100
96 : DUMPATOMS FILE=dump.xyz ATOMS=rna,ions STRIDE=100
97 : ```
98 :
99 : Notice that if the biased variable requires a molecule to be whole, you might have to put
100 : the [WHOLEMOLECULES](WHOLEMOLECULES.md) command before computing that variable and leave the default STRIDE=1.
101 :
102 : This command instructs plumed to center all atoms around the center of mass of a solute molecule.
103 :
104 : ```plumed
105 : solute: GROUP ATOMS=1-100
106 : all: GROUP ATOMS=1-1000
107 : # center of the solute:
108 : # notice that since plumed 2.2 this also works if the
109 : # solute molecule is broken
110 : com: COM ATOMS=solute
111 : # notice that we wrap around a single atom. this should be fast
112 : WRAPAROUND ATOMS=all AROUND=com
113 : DUMPATOMS FILE=dump.xyz ATOMS=all
114 : ```
115 :
116 : Notice that whereas [WHOLEMOLECULES](WHOLEMOLECULES.md) is designed to make molecules whole,
117 : WRAPAROUND can easily break molecules. In the last example,
118 : if solvent (atoms 101-1000) is made e.g. of water, then water
119 : molecules could be broken by WRAPAROUND (hydrogen could end up
120 : in an image and oxygen in another one).
121 : One solution is to use [WHOLEMOLECULES](WHOLEMOLECULES.md) on _all_ the water molecules
122 : after WRAPAROUND. This is tedious. A better solution is to use the
123 : GROUPBY option which is going
124 : to consider the atoms listed in ATOMS as a list of groups
125 : each of size GROUPBY. The first atom of the group will be brought
126 : close to the AROUND atoms. The following atoms of the group
127 : will be just brought close to the first atom of the group.
128 : Assuming that oxygen is the first atom of each water molecules,
129 : in the following examples all the water oxygen atoms will be brought
130 : close to the solute, and all the hydrogen atoms will be kept close
131 : to their related oxygen.
132 :
133 : ```plumed
134 : solute: GROUP ATOMS=1-100
135 : water: GROUP ATOMS=101-1000
136 : com: COM ATOMS=solute
137 : # notice that we wrap around a single atom. this should be fast
138 : WRAPAROUND ATOMS=solute AROUND=com
139 : # notice that we wrap around a single atom. this should be fast
140 : WRAPAROUND ATOMS=water AROUND=com GROUPBY=3
141 : DUMPATOMS FILE=dump.xyz ATOMS=solute,water
142 : ```
143 :
144 : */
145 : //+ENDPLUMEDOC
146 :
147 :
148 : class WrapAround:
149 : public ActionPilot,
150 : public ActionAtomistic {
151 : // cppcheck-suppress duplInheritedMember
152 : std::vector<Vector> refatoms;
153 : std::vector<std::pair<std::size_t,std::size_t> > p_atoms;
154 : std::vector<std::pair<std::size_t,std::size_t> > p_reference;
155 : unsigned groupby;
156 : bool pair_;
157 : public:
158 : explicit WrapAround(const ActionOptions&ao);
159 : static void registerKeywords( Keywords& keys );
160 0 : bool actionHasForces() override {
161 0 : return false;
162 : }
163 : void calculate() override;
164 590 : void apply() override {}
165 : };
166 :
167 : PLUMED_REGISTER_ACTION(WrapAround,"WRAPAROUND")
168 :
169 8 : void WrapAround::registerKeywords( Keywords& keys ) {
170 8 : Action::registerKeywords( keys );
171 8 : ActionAtomistic::registerKeywords( keys );
172 8 : ActionPilot::registerKeywords( keys );
173 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!");
174 8 : keys.add("atoms","AROUND","reference atoms");
175 8 : keys.add("atoms","ATOMS","wrapped atoms");
176 8 : keys.add("compulsory","GROUPBY","1","group atoms so as not to break molecules");
177 8 : keys.addFlag("PAIR", false, "Pair atoms in AROUND and ATOMS groups");
178 8 : }
179 :
180 6 : WrapAround::WrapAround(const ActionOptions&ao):
181 : Action(ao),
182 : ActionPilot(ao),
183 : ActionAtomistic(ao),
184 6 : groupby(1),
185 6 : pair_(false) {
186 : std::vector<AtomNumber> atoms;
187 12 : parseAtomList("ATOMS",atoms);
188 : std::vector<AtomNumber> reference;
189 6 : parseAtomList("AROUND",reference);
190 6 : parse("GROUPBY",groupby);
191 6 : parseFlag("PAIR", pair_);
192 :
193 6 : log.printf(" atoms in reference :");
194 13 : for(unsigned j=0; j<reference.size(); ++j) {
195 7 : log.printf(" %d",reference[j].serial() );
196 : }
197 6 : log.printf("\n");
198 6 : log.printf(" atoms to be wrapped :");
199 422 : for(unsigned j=0; j<atoms.size(); ++j) {
200 416 : log.printf(" %d",atoms[j].serial() );
201 : }
202 6 : log.printf("\n");
203 6 : if(groupby>1) {
204 1 : log<<" atoms will be grouped by "<<groupby<<"\n";
205 : }
206 6 : if(pair_) {
207 0 : log.printf(" pairing atoms and references\n");
208 : }
209 :
210 6 : if(atoms.size()%groupby!=0) {
211 0 : error("number of atoms should be a multiple of groupby option");
212 : }
213 : // additional checks with PAIR
214 6 : if(pair_ && atoms.size()!=reference.size()*groupby) {
215 0 : error("with PAIR you must have: #ATOMS = #AROUND * #GROUPBY");
216 : }
217 :
218 6 : checkRead();
219 :
220 : // do not remove duplicates with pair
221 6 : if(!pair_) {
222 6 : if(groupby<=1) {
223 5 : Tools::removeDuplicates(atoms);
224 : }
225 6 : Tools::removeDuplicates(reference);
226 : }
227 :
228 6 : std::vector<AtomNumber> merged(atoms.size()+reference.size());
229 6 : merge(atoms.begin(),atoms.end(),reference.begin(),reference.end(),merged.begin());
230 6 : p_atoms.resize( atoms.size() );
231 422 : for(unsigned i=0; i<atoms.size(); ++i) {
232 416 : p_atoms[i] = getValueIndices( atoms[i] );
233 : }
234 6 : refatoms.resize( reference.size() );
235 6 : p_reference.resize( reference.size() );
236 13 : for(unsigned i=0; i<reference.size(); ++i) {
237 7 : p_reference[i] = getValueIndices( reference[i] );
238 : }
239 6 : Tools::removeDuplicates(merged);
240 6 : requestAtoms(merged);
241 : doNotRetrieve();
242 : doNotForce();
243 6 : }
244 :
245 590 : void WrapAround::calculate() {
246 1185 : for(unsigned j=0; j<p_reference.size(); ++j) {
247 595 : refatoms[j] = getGlobalPosition(p_reference[j]);
248 : }
249 :
250 15265 : for(unsigned i=0; i<p_atoms.size(); i+=groupby) {
251 14675 : Vector second, first=getGlobalPosition(p_atoms[i]);
252 : double mindist2=std::numeric_limits<double>::max();
253 : int closest=-1;
254 14675 : if(pair_) {
255 0 : closest = i/groupby;
256 : } else {
257 29900 : for(unsigned j=0; j<p_reference.size(); ++j) {
258 15225 : second=refatoms[j];
259 : const Vector distance=pbcDistance(first,second);
260 15225 : const double distance2=modulo2(distance);
261 15225 : if(distance2<mindist2) {
262 : mindist2=distance2;
263 14955 : closest=j;
264 : }
265 : }
266 14675 : plumed_massert(closest>=0,"closest not found");
267 : }
268 14675 : second=refatoms[closest];
269 : // place first atom of the group
270 14675 : first=second+pbcDistance(second,first);
271 14675 : setGlobalPosition(p_atoms[i],first);
272 : // then place other atoms close to the first of the group
273 15170 : for(unsigned j=1; j<groupby; j++) {
274 495 : second=getGlobalPosition(p_atoms[i+j]);
275 495 : setGlobalPosition( p_atoms[i+j], first+pbcDistance(first,second) );
276 : }
277 : }
278 590 : }
279 :
280 :
281 :
282 : }
283 :
284 : }
|
