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