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 : It is similar to the ALIGN_ATOMS keyword of plumed1, and is needed since some
45 : MD dynamics code (e.g. GROMACS) can break molecules during the calculation.
46 :
47 : Running some CVs without this command can cause there to be discontinuities changes
48 : in the CV value and artifacts in the calculations. This command can be applied
49 : more than once. To see what effect is has use a variable without pbc or use
50 : the \ref DUMPATOMS directive to output the atomic positions.
51 :
52 : \attention
53 : This directive modifies the stored position at the precise moment
54 : it is executed. This means that only collective variables
55 : which are below it in the input script will see the corrected positions.
56 : As a general rule, put it at the top of the input file. Also, unless you
57 : know exactly what you are doing, leave the default stride (1), so that
58 : this action is performed at every MD step.
59 :
60 : The way WHOLEMOLECULES modifies each of the listed entities is this:
61 : - First atom of the list is left in place
62 : - Each atom of the list is shifted by a lattice vectors so that it becomes as close as possible
63 : to the previous one, iteratively.
64 :
65 : In this way, if an entity consists of a list of atoms such that consecutive atoms in the
66 : list are always closer than half a box side the entity will become whole.
67 : This can be usually achieved selecting consecutive atoms (1-100), but it is also possible
68 : to skip some atoms, provided consecutive chosen atoms are close enough.
69 :
70 : \par Examples
71 :
72 : This command instructs plumed to reconstruct the molecule containing atoms 1-20
73 : at every step of the calculation and dump them on a file.
74 :
75 : \plumedfile
76 : # to see the effect, one could dump the atoms as they were before molecule reconstruction:
77 : # DUMPATOMS FILE=dump-broken.xyz ATOMS=1-20
78 : WHOLEMOLECULES ENTITY0=1-20
79 : DUMPATOMS FILE=dump.xyz ATOMS=1-20
80 : \endplumedfile
81 :
82 : This command instructs plumed to reconstruct two molecules containing atoms 1-20 and 30-40
83 :
84 : \plumedfile
85 : WHOLEMOLECULES ENTITY0=1-20 ENTITY1=30-40
86 : DUMPATOMS FILE=dump.xyz ATOMS=1-20,30-40
87 : \endplumedfile
88 :
89 : This command instructs plumed to reconstruct the chain of backbone atoms in a
90 : protein
91 :
92 : \plumedfile
93 : #SETTINGS MOLFILE=regtest/basic/rt32/helix.pdb
94 : MOLINFO STRUCTURE=helix.pdb
95 : WHOLEMOLECULES RESIDUES=all MOLTYPE=protein
96 : \endplumedfile
97 :
98 : */
99 : //+ENDPLUMEDOC
100 :
101 :
102 : class WholeMolecules:
103 : public ActionPilot,
104 : public ActionAtomistic {
105 : std::vector<std::vector<std::pair<std::size_t,std::size_t> > > p_groups;
106 : std::vector<std::vector<std::pair<std::size_t,std::size_t> > > p_roots;
107 : std::vector<Vector> refs;
108 : bool doemst, addref;
109 : public:
110 : explicit WholeMolecules(const ActionOptions&ao);
111 : static void registerKeywords( Keywords& keys );
112 841 : bool actionHasForces() override {
113 841 : return false;
114 : }
115 : void calculate() override;
116 11946 : void apply() override {}
117 : };
118 :
119 : PLUMED_REGISTER_ACTION(WholeMolecules,"WHOLEMOLECULES")
120 :
121 95 : void WholeMolecules::registerKeywords( Keywords& keys ) {
122 95 : Action::registerKeywords( keys );
123 95 : ActionPilot::registerKeywords( keys );
124 95 : ActionAtomistic::registerKeywords( keys );
125 190 : 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!");
126 190 : 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,...");
127 190 : keys.reset_style("ENTITY","atoms");
128 190 : 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 "
129 : "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 "
130 : "specifying all. Alternatively, if you wish to use a subset of the residues you can specify the particular residues "
131 : "you are interested in as a list of numbers");
132 190 : keys.add("optional","MOLTYPE","the type of molecule that is under study. This is used to define the backbone atoms");
133 190 : keys.addFlag("EMST", false, "Define atoms sequence in entities using an Euclidean minimum spanning tree");
134 190 : keys.addFlag("ADDREFERENCE", false, "Define the reference position of the first atom of each entity using a PDB file");
135 95 : }
136 :
137 71 : WholeMolecules::WholeMolecules(const ActionOptions&ao):
138 : Action(ao),
139 : ActionPilot(ao),
140 : ActionAtomistic(ao),
141 71 : doemst(false), addref(false) {
142 : std::vector<std::vector<AtomNumber> > groups;
143 : std::vector<std::vector<AtomNumber> > roots;
144 : // parse optional flags
145 71 : parseFlag("EMST", doemst);
146 142 : parseFlag("ADDREFERENCE", addref);
147 :
148 : // create groups from ENTITY
149 422 : for(int i=0;; i++) {
150 : std::vector<AtomNumber> group;
151 986 : parseAtomList("ENTITY",i,group);
152 493 : if( group.empty() ) {
153 : break;
154 : }
155 422 : groups.push_back(group);
156 422 : }
157 :
158 : // Read residues to align from MOLINFO
159 : std::vector<std::string> resstrings;
160 142 : parseVector("RESIDUES",resstrings);
161 71 : if( resstrings.size()>0 ) {
162 0 : if( resstrings.size()==1 ) {
163 0 : if( resstrings[0]=="all" ) {
164 : resstrings[0]="all-ter"; // Include terminal groups in alignment
165 : }
166 : }
167 : std::string moltype;
168 0 : parse("MOLTYPE",moltype);
169 0 : if(moltype.length()==0) {
170 0 : error("Found RESIDUES keyword without specification of the molecule - use MOLTYPE");
171 : }
172 0 : auto* moldat=plumed.getActionSet().selectLatest<GenericMolInfo*>(this);
173 0 : if( !moldat ) {
174 0 : error("MOLINFO is required to use RESIDUES");
175 : }
176 : std::vector< std::vector<AtomNumber> > backatoms;
177 0 : moldat->getBackbone( resstrings, moltype, backatoms );
178 0 : for(unsigned i=0; i<backatoms.size(); ++i) {
179 0 : groups.push_back( backatoms[i] );
180 : }
181 0 : }
182 :
183 : // check number of groups
184 71 : if(groups.size()==0) {
185 0 : error("no atoms found for WHOLEMOLECULES!");
186 : }
187 :
188 : // if using PDBs reorder atoms in groups based on proximity in PDB file
189 71 : if(doemst) {
190 2 : auto* moldat=plumed.getActionSet().selectLatest<GenericMolInfo*>(this);
191 2 : if( !moldat ) {
192 0 : error("MOLINFO is required to use EMST");
193 : }
194 : // initialize tree
195 2 : Tree tree = Tree(moldat);
196 : // cycle on groups and reorder atoms
197 4 : for(unsigned i=0; i<groups.size(); ++i) {
198 4 : groups[i] = tree.getTree(groups[i]);
199 : // store root atoms
200 4 : roots.push_back(tree.getRoot());
201 : }
202 : } else {
203 : // fill root vector with previous atom in groups
204 489 : for(unsigned i=0; i<groups.size(); ++i) {
205 : std::vector<AtomNumber> root;
206 11390 : for(unsigned j=0; j<groups[i].size()-1; ++j) {
207 10970 : root.push_back(groups[i][j]);
208 : }
209 : // store root atoms
210 420 : roots.push_back(root);
211 : }
212 : }
213 :
214 : // adding reference if needed
215 71 : if(addref) {
216 2 : auto* moldat=plumed.getActionSet().selectLatest<GenericMolInfo*>(this);
217 2 : if( !moldat ) {
218 0 : error("MOLINFO is required to use ADDREFERENCE");
219 : }
220 4 : for(unsigned i=0; i<groups.size(); ++i) {
221 : // add reference position of first atom in entity
222 4 : refs.push_back(moldat->getPosition(groups[i][0]));
223 : }
224 : }
225 :
226 : // print out info
227 493 : for(unsigned i=0; i<groups.size(); ++i) {
228 422 : log.printf(" atoms in entity %d : ",i);
229 12057 : for(unsigned j=0; j<groups[i].size(); ++j) {
230 11635 : log.printf("%d ",groups[i][j].serial() );
231 : }
232 422 : log.printf("\n");
233 422 : if(addref) {
234 2 : log.printf(" with reference position : %lf %lf %lf\n",refs[i][0],refs[i][1],refs[i][2]);
235 : }
236 : }
237 :
238 : // collect all atoms
239 : std::vector<AtomNumber> merge;
240 493 : for(unsigned i=0; i<groups.size(); ++i) {
241 422 : merge.insert(merge.end(),groups[i].begin(),groups[i].end());
242 : }
243 :
244 : // Convert groups to p_groups
245 71 : p_groups.resize( groups.size() );
246 493 : for(unsigned i=0; i<groups.size(); ++i) {
247 422 : p_groups[i].resize( groups[i].size() );
248 12057 : for(unsigned j=0; j<groups[i].size(); ++j) {
249 11635 : p_groups[i][j] = getValueIndices( groups[i][j] );
250 : }
251 : }
252 : // Convert roots to p_roots
253 71 : p_roots.resize( roots.size() );
254 493 : for(unsigned i=0; i<roots.size(); ++i) {
255 422 : p_roots[i].resize( roots[i].size() );
256 11635 : for(unsigned j=0; j<roots[i].size(); ++j) {
257 11213 : p_roots[i][j] = getValueIndices( roots[i][j] );
258 : }
259 : }
260 :
261 :
262 71 : checkRead();
263 71 : Tools::removeDuplicates(merge);
264 71 : requestAtoms(merge);
265 : doNotRetrieve();
266 : doNotForce();
267 71 : }
268 :
269 11946 : void WholeMolecules::calculate() {
270 24596 : for(unsigned i=0; i<p_groups.size(); ++i) {
271 12650 : Vector first = getGlobalPosition(p_groups[i][0]);
272 12650 : if(addref) {
273 12 : first = refs[i]+pbcDistance(refs[i],first);
274 12 : setGlobalPosition( p_groups[i][0], first );
275 : }
276 12650 : if(!doemst) {
277 662472 : for(unsigned j=1; j<p_groups[i].size(); ++j) {
278 649826 : Vector second=getGlobalPosition(p_groups[i][j]);
279 649826 : first = first+pbcDistance(first,second);
280 649826 : setGlobalPosition(p_groups[i][j], first );
281 : }
282 : } else {
283 490 : for(unsigned j=1; j<p_groups[i].size(); ++j) {
284 486 : Vector first=getGlobalPosition(p_roots[i][j-1]);
285 486 : Vector second=getGlobalPosition(p_groups[i][j]);
286 486 : second=first+pbcDistance(first,second);
287 486 : setGlobalPosition(p_groups[i][j], second );
288 : }
289 : }
290 : }
291 11946 : }
292 :
293 :
294 : }
295 : }
|
