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