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