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