Line data Source code
1 : /* +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
2 : Copyright (c) 2011-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 : #include "tools/OpenMP.h"
32 : #include "tools/Tree.h"
33 :
34 : #include <vector>
35 : #include <string>
36 :
37 : namespace PLMD {
38 : namespace generic {
39 :
40 : //+PLUMEDOC GENERIC WHOLEMOLECULES
41 : /*
42 : This action is used to rebuild molecules that can become split by the periodic boundary conditions.
43 :
44 : This command performs an operation that is similar what was done by the ALIGN_ATOMS keyword from plumed1.
45 : This operation is needed as some MD dynamics code (e.g. GROMACS) can break molecules during the calculation.
46 : Whenever we are able we try to ensure that molecules are reconstructed automatically. You thus do not need
47 : to use this action when you are using actions such as [COM](COM.md), [CENTER](CENTER.md), [GYRATION](GYRATION.md)
48 : and so on as the reconstruction of molecules is done automatically in these actions. It is, however, important to understand
49 : molecule reconstruction as there are many cases (e.g. when you are calculating the end-to-end distance of a polymer)
50 : where not using the WHOLEMOLCULES command can cause the artifacts discussed in the attached reference.
51 :
52 : If you think that you need to use this command a good idea is to use the [DUMPATOMS](DUMPATOMS.md) directive
53 : to output the atomic positions. This will allow you to see the effect that including/not including WHOLEMOLECULES
54 : has on the calculation.
55 :
56 : !!! attention "modifies stored positions"
57 :
58 : This directive modifies the stored position at the precise moment
59 : it is executed. This means that only collective variables
60 : which are below it in the input script will see the corrected positions.
61 : As a general rule, put it at the top of the input file. Also, unless you
62 : know exactly what you are doing, leave the default stride (1), so that
63 : this action is performed at every MD step.
64 :
65 : Notice that the behavior of WHOLEMOLECULES is affected by the last [MOLINFO](MOLINFO.md) action
66 : that is present in the input file before the WHOLEMOLECULES command. Specifically, if the
67 : [MOLINFO](MOLINFO.md) action does not have a `WHOLE` flag, then the behavior is the following:
68 :
69 : - The first atom of the list is left in place
70 : - Each atom of the list is shifted by a lattice vectors so that it becomes as close as possible
71 : to the previous one, iteratively.
72 :
73 : In this way, if an entity consists of a list of atoms such that consecutive atoms in the
74 : list are always closer than half a box side the entity will become whole.
75 : This can be usually achieved selecting consecutive atoms (1-100), but it is also possible
76 : to skip some atoms, provided consecutive chosen atoms are close enough.
77 :
78 : If, by contrast, the [MOLINFO](MOLINFO.md) action does have a `WHOLE` flag, then a minimum spanning tree
79 : is built based on the atoms passed to WHOLEMOLECULES using the coordinates in the PDB
80 : passed to [MOLINFO](MOLINFO.md) as a reference, and this tree is used to reconstruct PBCs.
81 : This approach is more robust when dealing with complexes of multiple molecules.
82 :
83 : ## Examples
84 :
85 : This command instructs plumed to reconstruct the molecule containing atoms 1-20
86 : at every step of the calculation and dump them on a file.
87 :
88 : ```plumed
89 : # to see the effect, one could dump the atoms as they were before molecule reconstruction:
90 : # DUMPATOMS FILE=dump-broken.xyz ATOMS=1-20
91 : WHOLEMOLECULES ENTITY0=1-20
92 : DUMPATOMS FILE=dump.xyz ATOMS=1-20
93 : ```
94 :
95 : This command instructs plumed to reconstruct two molecules containing atoms 1-20 and 30-40
96 :
97 : ```plumed
98 : WHOLEMOLECULES ENTITY0=1-20 ENTITY1=30-40
99 : DUMPATOMS FILE=dump.xyz ATOMS=1-20,30-40
100 : ```
101 :
102 : This command instructs plumed to reconstruct the chain of backbone atoms in a
103 : protein
104 :
105 : ```plumed
106 : #SETTINGS MOLFILE=regtest/basic/rt32/helix.pdb
107 : MOLINFO STRUCTURE=regtest/basic/rt32/helix.pdb
108 : WHOLEMOLECULES RESIDUES=all MOLTYPE=protein
109 : ```
110 :
111 : */
112 : //+ENDPLUMEDOC
113 :
114 :
115 : class WholeMolecules:
116 : public ActionPilot,
117 : public ActionAtomistic {
118 : std::vector<std::vector<std::pair<std::size_t,std::size_t> > > p_groups;
119 : std::vector<std::vector<std::pair<std::size_t,std::size_t> > > p_roots;
120 : std::vector<Vector> refs= {};
121 : bool doemst=false;
122 : bool addref=false;
123 : public:
124 : explicit WholeMolecules(const ActionOptions&ao);
125 : static void registerKeywords( Keywords& keys );
126 841 : bool actionHasForces() override {
127 841 : return false;
128 : }
129 : void calculate() override;
130 13668 : void apply() override {}
131 : };
132 :
133 : PLUMED_REGISTER_ACTION(WholeMolecules,"WHOLEMOLECULES")
134 :
135 99 : void WholeMolecules::registerKeywords( Keywords& keys ) {
136 99 : Action::registerKeywords( keys );
137 99 : ActionPilot::registerKeywords( keys );
138 99 : ActionAtomistic::registerKeywords( keys );
139 99 : 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!");
140 99 : keys.add("numbered","ENTITY","the atoms that make up a molecule that you wish to align. To specify multiple molecules use a list of ENTITY keywords: ENTITY0, ENTITY1,...");
141 198 : keys.reset_style("ENTITY","atoms");
142 99 : keys.add("residues","RESIDUES","this command specifies that the backbone atoms in a set of residues all must be aligned. It must be used in tandem with the \\ref MOLINFO "
143 : "action and the MOLTYPE keyword. If you wish to use all the residues from all the chains in your system you can do so by "
144 : "specifying all. Alternatively, if you wish to use a subset of the residues you can specify the particular residues "
145 : "you are interested in as a list of numbers");
146 99 : keys.add("optional","MOLTYPE","the type of molecule that is under study. This is used to define the backbone atoms");
147 99 : keys.addFlag("EMST", false, "only for backward compatibility, as of PLUMED 2.11 this is the default when using MOLINFO with WHOLE");
148 99 : keys.addFlag("ADDREFERENCE", false, "Define the reference position of the first atom of each entity using a PDB file");
149 99 : keys.addDOI("10.1007/978-1-4939-9608-7_21");
150 99 : }
151 :
152 76 : WholeMolecules::WholeMolecules(const ActionOptions&ao):
153 : Action(ao),
154 : ActionPilot(ao),
155 76 : ActionAtomistic(ao) {
156 : std::vector<std::vector<AtomNumber> > groups;
157 : std::vector<std::vector<AtomNumber> > roots;
158 : // parse optional flags
159 : bool doemst_tmp;
160 76 : parseFlag("EMST", doemst_tmp);
161 76 : if(doemst_tmp) {
162 1 : log << "EMST option is not needed any more as of PLUMED 2.11\n";
163 : }
164 76 : parseFlag("ADDREFERENCE", addref);
165 :
166 76 : auto* infomoldat=plumed.getActionSet().selectLatest<GenericMolInfo*>(this);
167 :
168 : // create groups from ENTITY
169 433 : for(int i=0;; i++) {
170 : std::vector<AtomNumber> group;
171 1018 : parseAtomList("ENTITY",i,group);
172 509 : if( group.empty() ) {
173 : break;
174 : }
175 433 : groups.push_back(group);
176 433 : }
177 :
178 : // Read residues to align from MOLINFO
179 : std::vector<std::string> resstrings;
180 152 : parseVector("RESIDUES",resstrings);
181 76 : if( resstrings.size()>0 ) {
182 0 : if( resstrings.size()==1 ) {
183 0 : if( resstrings[0]=="all" ) {
184 : resstrings[0]="all-ter"; // Include terminal groups in alignment
185 : }
186 : }
187 : std::string moltype;
188 0 : parse("MOLTYPE",moltype);
189 0 : if(moltype.length()==0) {
190 0 : error("Found RESIDUES keyword without specification of the molecule - use MOLTYPE");
191 : }
192 0 : if( !infomoldat ) {
193 0 : error("MOLINFO is required to use RESIDUES");
194 : }
195 : std::vector< std::vector<AtomNumber> > backatoms;
196 0 : infomoldat->getBackbone( resstrings, moltype, backatoms );
197 0 : for(unsigned i=0; i<backatoms.size(); ++i) {
198 0 : groups.push_back( backatoms[i] );
199 : }
200 0 : }
201 :
202 : // check number of groups
203 76 : if(groups.size()==0) {
204 0 : error("no atoms found for WHOLEMOLECULES!");
205 : }
206 :
207 : // if using PDBs reorder atoms in groups based on proximity in PDB file
208 76 : if(infomoldat && infomoldat->isWhole()) {
209 2 : doemst=true;
210 : }
211 :
212 76 : if(doemst_tmp && ! doemst) {
213 0 : error("cannot enable EMST if MOLINFO is not WHOLE");
214 : }
215 :
216 76 : if(doemst) {
217 2 : if( !infomoldat ) {
218 0 : error("MOLINFO is required to use EMST");
219 : }
220 : // initialize tree
221 2 : Tree myTree = Tree(infomoldat);
222 : // cycle on groups and reorder atoms
223 4 : for(unsigned i=0; i<groups.size(); ++i) {
224 2 : groups[i] = myTree.getTree(groups[i]);
225 : // store root atoms
226 2 : roots.push_back(myTree.getRoot());
227 : }
228 2 : } else {
229 : // fill root vector with previous atom in groups
230 505 : for(unsigned i=0; i<groups.size(); ++i) {
231 : std::vector<AtomNumber> root;
232 11885 : for(unsigned j=0; j<groups[i].size()-1; ++j) {
233 11454 : root.push_back(groups[i][j]);
234 : }
235 : // store root atoms
236 431 : roots.push_back(root);
237 : }
238 : }
239 :
240 : // adding reference if needed
241 76 : if(addref) {
242 2 : if( !infomoldat ) {
243 0 : error("MOLINFO is required to use ADDREFERENCE");
244 : }
245 4 : for(unsigned i=0; i<groups.size(); ++i) {
246 : // add reference position of first atom in entity
247 4 : refs.push_back(infomoldat->getPosition(groups[i][0]));
248 : }
249 : }
250 :
251 : // print out info
252 509 : for(unsigned i=0; i<groups.size(); ++i) {
253 433 : log.printf(" atoms in entity %d : ",i);
254 12563 : for(unsigned j=0; j<groups[i].size(); ++j) {
255 12130 : log.printf("%d ",groups[i][j].serial() );
256 : }
257 433 : log.printf("\n");
258 433 : if(addref) {
259 2 : log.printf(" with reference position : %lf %lf %lf\n",refs[i][0],refs[i][1],refs[i][2]);
260 : }
261 : }
262 :
263 : // collect all atoms
264 : std::vector<AtomNumber> merge;
265 509 : for(unsigned i=0; i<groups.size(); ++i) {
266 433 : merge.insert(merge.end(),groups[i].begin(),groups[i].end());
267 : }
268 :
269 : // Convert groups to p_groups
270 76 : p_groups.resize( groups.size() );
271 509 : for(unsigned i=0; i<groups.size(); ++i) {
272 433 : p_groups[i].resize( groups[i].size() );
273 12563 : for(unsigned j=0; j<groups[i].size(); ++j) {
274 12130 : p_groups[i][j] = getValueIndices( groups[i][j] );
275 : }
276 : }
277 : // Convert roots to p_roots
278 76 : p_roots.resize( roots.size() );
279 509 : for(unsigned i=0; i<roots.size(); ++i) {
280 433 : p_roots[i].resize( roots[i].size() );
281 12130 : for(unsigned j=0; j<roots[i].size(); ++j) {
282 11697 : p_roots[i][j] = getValueIndices( roots[i][j] );
283 : }
284 : }
285 :
286 76 : checkRead();
287 76 : Tools::removeDuplicates(merge);
288 76 : requestAtoms(merge);
289 : doNotRetrieve();
290 : doNotForce();
291 76 : }
292 :
293 13668 : void WholeMolecules::calculate() {
294 31562 : for(unsigned i=0; i<p_groups.size(); ++i) {
295 17894 : Vector first = getGlobalPosition(p_groups[i][0]);
296 17894 : if(addref) {
297 12 : first = refs[i]+pbcDistance(refs[i],first);
298 12 : setGlobalPosition( p_groups[i][0], first );
299 : }
300 17894 : if(!doemst) {
301 896248 : for(unsigned j=1; j<p_groups[i].size(); ++j) {
302 878358 : Vector second=getGlobalPosition(p_groups[i][j]);
303 878358 : first = first+pbcDistance(first,second);
304 878358 : setGlobalPosition(p_groups[i][j], first );
305 : }
306 : } else {
307 490 : for(unsigned j=1; j<p_groups[i].size(); ++j) {
308 486 : Vector firstPos=getGlobalPosition(p_roots[i][j-1]);
309 486 : Vector secondPos=getGlobalPosition(p_groups[i][j]);
310 486 : secondPos=firstPos+pbcDistance(firstPos,secondPos);
311 486 : setGlobalPosition(p_groups[i][j], secondPos );
312 : }
313 : }
314 : }
315 13668 : }
316 :
317 :
318 : }
319 : }
|
