Line data Source code
1 : /* +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
2 : Copyright (c) 2019-2020
3 :
4 : This file is part of funnel code module.
5 :
6 : The FM code respects the CC BY-NC license.
7 : Users are free to download, adapt and use the code as long as it is not for commercial purposes.
8 : +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ */
9 : #include "colvar/Colvar.h"
10 : #include "core/ActionRegister.h"
11 : #include <string>
12 : #include <cmath>
13 : #include <cassert>
14 : #include <vector>
15 : #include "tools/Tools.h"
16 : #include "tools/PDB.h"
17 : #include <iostream>
18 : #include "tools/RMSD.h"
19 :
20 : using namespace std;
21 :
22 : namespace PLMD {
23 : namespace funnel {
24 :
25 : //+PLUMEDOC COLVAR FUNNEL_PS
26 : /*
27 : FUNNEL_PS implements the Funnel-Metadynamics (FM) technique in PLUMED 2.
28 :
29 : Please read the first two papers in the references below to better understand the notions hereby reported.
30 :
31 : This colvar evaluates the position of a ligand of interest with respect to a given line, built from the two points
32 : A and B, and should be used together with the [FUNNEL](FUNNEL.md) bias.
33 : The constructed line represents the axis of the funnel-shape restraint potential, which should be placed so
34 : as to include the portion of a macromolecule (i.e., protein, DNA, etc.) that should be explored.
35 : Since it is important that the position of the line is updated based on the motion of the macromolecule during
36 : the simulation, this colvar incorporates an alignment method. The latter uses the TYPE=OPTIMAL option to remove
37 : motions due to rotation and translation of the macromolecule with respect to a reference structure, which is
38 : provided by the user. In order to accomplish the task, an optimal alignment matrix is calculated using the
39 : Kearsley algorithm that is discussed in the third paper cited below.
40 : The reference structure should be given as a pdb file, containing only the atoms of the macromolecule or a
41 : selection of them (e.g., the protein CA atoms of secondary structures). In contrast to the methods discussed in
42 : the documentation for [RMSD](RMSD.md), the values reported in the occupancy and beta-factor columns of the pdb file are not important
43 : since they will be overwritten and replaced by the value 1.00 during the procedure. It is important to understand
44 : that all atoms in the file will be used for the alignment, even if they display 0.00 in the occupancy column.
45 :
46 : The ligand can be represented by one single atom or the center of mass (COM) of a group of atoms that should be
47 : provided by the user.
48 :
49 : By default FUNNEL_PS is computed taking into account periodic boundary conditions. Since PLUMED 2.5, molecules are
50 : rebuilt using a procedure that is equivalent to that done in [WHOLEMOLECULES](WHOLEMOLECULES.md). We note that this action is local
51 : to this colvar, thus it does not modify the coordinates stored in PLUMED. Moreover, FUNNEL_PS requires an ANCHOR atom
52 : to be specified in order to facilitate the reconstruction of periodic boundary conditions. This atom should be the
53 : closest macromolecule's atom to the ligand and it should reduce the risk of ligand "warping" in the simulation box.
54 : Nevertheless, we highly recommend to add to the PLUMED input file a custom line of [WHOLEMOLECULES](WHOLEMOLECULES.md), in order to
55 : be sure of reconstructing the ligand together with the macromolecule (please look the examples). In this case, the user
56 : can use the NOPBC flag to turn off the internal periodic boundary condition reconstruction.
57 :
58 : FUNNEL_PS is divided in two components (fps.lp and fps.ld) which evaluate the projection of the ligand along the funnel line
59 : and the distance from it, respectively. The values attributed to these two components are then used together with the
60 : potential file created by the [FUNNEL](FUNNEL.md) bias to define if the ligand is within or not in the funnel-shape restraint
61 : potential. In the latter case, the potential will force the ligand to enter within the funnel boundaries.
62 :
63 : ## Examples
64 :
65 : The following input tells plumed to print the FUNNEL_PS components for the COM of a ligand. The inputs are a reference structure,
66 : which is the structure used for the alignment, the COM of the molecule we want to track, and 2 points in the Cartesian space
67 : (i.e., x, y, and z) to draw the line representing the funnel axis.
68 :
69 : ```plumed
70 : lig: COM ATOMS=2446,2447,2448,2449,2451
71 : fps: FUNNEL_PS REFERENCE=protein.pdb LIGAND=lig POINTS=5.3478,-0.7278,2.4746,7.3785,6.7364,-9.3624
72 : PRINT ARG=fps.lp,fps.ld
73 : ```
74 :
75 : It is recommended to add a line to force the reconstruction of the periodic boundary conditions. In the following example,
76 : [WHOLEMOLECULES](WHOLEMOLECULES.md) was added to make sure that the ligand was reconstructed together with the protein. The list contains
77 : all the atoms reported in the start.pdb file followed by the ANCHOR atom and the ligand. All atoms should be contained in the
78 : same entity and the correct order.
79 :
80 : ```plumed
81 : WHOLEMOLECULES ENTITY0=54,75,212,228,239,258,311,328,348,372,383,402,421,463,487,503,519,657,674,690,714,897,914,934,953,964,
82 : 974,985,1007,1018,1037,1247,1264,1283,1302,1324,1689,1708,1727,1738,1962,1984,1994,2312,2329,2349,2467,2490,2500,2517,2524,2536,
83 : 2547,2554,2569,2575,2591,2607,2635,2657,2676,2693,2700,2719,2735,2746,2770,2777,2788,2795,2805,2815,2832,2854,2868,2898,2904,
84 : 2911,2927,2948,2962,2472,3221,3224,3225,3228,3229,3231,3233,3235,3237
85 : lig: COM ATOMS=3221,3224,3225,3228,3229,3231,3233,3235,3237
86 : fps: FUNNEL_PS LIGAND=lig REFERENCE=start.pdb ANCHOR=2472 POINTS=4.724,5.369,4.069,4.597,5.721,4.343
87 : PRINT ARG=fps.lp,fps.ld
88 : ```
89 :
90 : */
91 : //+ENDPLUMEDOC
92 :
93 :
94 : class FUNNEL_PS : public Colvar {
95 :
96 : // Those are arrays that I created for the colvar
97 : std::vector<AtomNumber> ligand_com;
98 : std::vector<AtomNumber> anchor;
99 : std::vector<AtomNumber> numbers;
100 : bool pbc;
101 : PLMD::RMSD alignment;
102 : PLMD::PDB pdb;
103 : bool squared;
104 : private:
105 : vector<double> points;
106 : public:
107 : explicit FUNNEL_PS(const ActionOptions&);
108 : // active methods:
109 : virtual void calculate();
110 : static void registerKeywords(Keywords& keys);
111 : // I need a method in RMSDCoreCalc and these were requested
112 : std::vector<double> align;
113 : std::vector<double> displace;
114 : };
115 :
116 : using namespace std;
117 :
118 : PLUMED_REGISTER_ACTION(FUNNEL_PS,"FUNNEL_PS")
119 :
120 9 : void FUNNEL_PS::registerKeywords(Keywords& keys) {
121 9 : Colvar::registerKeywords( keys );
122 9 : keys.add("compulsory","REFERENCE","a file in pdb format containing the structure you would like to align.");
123 9 : keys.add("atoms","LIGAND","This MUST be a single atom, normally the COM of the ligand");
124 9 : keys.add("atoms","ANCHOR","Closest protein atom to the ligand, picked to avoid pbc problems during the simulation");
125 9 : keys.add("compulsory","POINTS","6 values defining x, y, and z of the 2 points used to construct the line. The order should be A_x,A_y,A_z,B_x,B_y,B_z.");
126 9 : keys.addFlag("SQUARED-ROOT",false,"Used to initialize the creation of the alignment variable");
127 18 : keys.addOutputComponent("lp","default","scalar","the position along the funnel line");
128 18 : keys.addOutputComponent("ld","default","scalar","the distance from the funnel line");
129 9 : keys.addDOI("10.1073/pnas.1303186110");
130 9 : keys.addDOI("10.1038/s41596-020-0342-4");
131 9 : keys.addDOI("10.1107/S0108767388010128");
132 9 : }
133 :
134 5 : FUNNEL_PS::FUNNEL_PS(const ActionOptions&ao):
135 : PLUMED_COLVAR_INIT(ao),
136 5 : pbc(true),squared(true) {
137 : string reference;
138 10 : parse("REFERENCE",reference);
139 : string type;
140 5 : type.assign("OPTIMAL");
141 10 : parseAtomList("LIGAND",ligand_com);
142 5 : if(ligand_com.size()!=1) {
143 0 : error("Number of specified atoms should be one, normally the COM of the ligand");
144 : }
145 5 : parseVector("POINTS",points);
146 5 : bool nopbc=!pbc;
147 5 : parseFlag("NOPBC",nopbc);
148 5 : pbc=!nopbc;
149 : bool sq;
150 5 : parseFlag("SQUARED-ROOT",sq);
151 5 : if (sq) {
152 0 : squared=false;
153 : }
154 5 : parseAtomList("ANCHOR",anchor);
155 5 : checkRead();
156 :
157 : // read everything in ang and transform to nm if we are not in natural units
158 5 : if( !pdb.read(reference,usingNaturalUnits(),0.1/getUnits().getLength()) ) {
159 0 : error("missing input file " + reference );
160 : }
161 :
162 : bool remove_com=true;
163 : bool normalize_weights=true;
164 : // here displace is a simple vector of ones
165 5 : align=pdb.getOccupancy();
166 16195 : for(unsigned i=0; i<align.size(); i++) {
167 16190 : align[i]=1.;
168 : } ;
169 5 : displace=pdb.getBeta();
170 16195 : for(unsigned i=0; i<displace.size(); i++) {
171 16190 : displace[i]=1.;
172 : } ;
173 : // reset again to reimpose uniform weights (safe to disable this)
174 5 : alignment.set(align,displace,pdb.getPositions(),type,remove_com,normalize_weights);
175 :
176 :
177 :
178 : // Array with inside both the structure to align and the atom to be aligned
179 5 : numbers=pdb.getAtomNumbers();
180 5 : numbers.push_back(anchor[0]);
181 5 : numbers.push_back(ligand_com[0]);
182 :
183 5 : log.printf(" average from file %s\n",reference.c_str());
184 5 : log.printf(" method for alignment : %s \n",type.c_str() );
185 :
186 5 : if(pbc) {
187 5 : log.printf(" using periodic boundary conditions\n");
188 : } else {
189 0 : log.printf(" without periodic boundary conditions\n");
190 : }
191 :
192 10 : addComponentWithDerivatives("lp");
193 10 : componentIsNotPeriodic("lp");
194 10 : addComponentWithDerivatives("ld");
195 5 : componentIsNotPeriodic("ld");
196 :
197 5 : requestAtoms( numbers );
198 :
199 5 : }
200 :
201 : // calculator
202 150 : void FUNNEL_PS::calculate() {
203 :
204 150 : if(pbc) {
205 150 : makeWhole();
206 : }
207 :
208 150 : Tensor Rotation;
209 150 : Matrix<std::vector<Vector> > drotdpos(3,3);
210 150 : std::vector<Vector> alignedpos;
211 150 : std::vector<Vector> centeredpos;
212 150 : std::vector<Vector> centeredref;
213 150 : std::vector<Vector> ddistdpos;
214 :
215 150 : std::vector<Vector> buffer;
216 :
217 : Vector centerreference;
218 150 : Vector centerpositions;
219 :
220 : // Created only to give the correct object to calc_FitElements
221 150 : std::vector<Vector> sourceAllPositions;
222 150 : std::vector<Vector> sourcePositions;
223 :
224 : // SourcePositions contains only the coordinates of the COM of the ligand, we need only this point
225 150 : sourceAllPositions=getPositions();
226 150 : sourcePositions=sourceAllPositions;
227 150 : sourcePositions.resize(sourcePositions.size()-2);// just the protein
228 :
229 : // The two points that define the axis : this can be moved in the initialization
230 150 : Vector p1 = VectorGeneric<3>(points[0],points[1],points[2]);
231 150 : Vector p2 = VectorGeneric<3>(points[3],points[4],points[5]);
232 : Vector s = p2 - p1;
233 :
234 : // I call the method calc_FitElements that initializes all feature that I need
235 : // except for centerreference that I need to calculate from scratch
236 : // Buffer has no meaning but I had to fulfill the requirements of calc_FitElements
237 150 : alignment.calc_FitElements( sourcePositions, Rotation, drotdpos, buffer, centerpositions, squared);
238 :
239 : // To Plumed developers: it would be interesting to make the functions to calculate centers of mass public or protected
240 : centerreference.zero();
241 485850 : for(unsigned i=0; i<pdb.size(); i++) {
242 485700 : centerreference+=pdb.getPositions()[i]*align[i]/align.size();
243 : }
244 :
245 : /*
246 : // I cancelled the additional lines in the library of RMSD.h, thus I am missing the center of the reference
247 : // Creating variable kito to extract only the center of the reference, since no method is calling
248 : // function getReferenceCenter()
249 : PLMD::RMSDCoreData* kito; kito = new RMSDCoreData(align,displace,sourcePositions,pdb->getPositions());
250 : centerreference = kito->getReferenceCenter();
251 : delete kito;
252 : */
253 :
254 : // DEBUG
255 :
256 : // the position is Rot(ligand-com_prot_md)+ com_prot_ref
257 : Vector ligand_centered =getPositions().back()-centerpositions;
258 150 : Vector ligand_aligned =matmul(Rotation,ligand_centered);
259 : Vector v = ligand_aligned +centerreference -p1;
260 :
261 : // DEBUG
262 : // log.printf(" ligando: %13.6lf %13.6lf %13.6lf\n",getPositions().back()[0],getPositions().back()[1],getPositions().back()[2] );
263 :
264 : //Projection vector v onto s
265 :
266 150 : Vector prj = (dotProduct(s,v)/dotProduct(s,s))*s;
267 :
268 : Vector height = v - prj;
269 150 : const double prj_height = height.modulo() ;
270 150 : const double inv_prj_height = 1.0/prj_height;
271 :
272 : // derivative of the prj: only on the com of the ligand
273 : Vector der_prj;
274 150 : der_prj=s/s.modulo();
275 :
276 : // derivative of the height: only on the com of the ligand
277 150 : Vector der_height(inv_prj_height*(height[0]-(s[0]/s.modulo2())*dotProduct(height,s)),
278 150 : inv_prj_height*(height[1]-(s[1]/s.modulo2())*dotProduct(height,s)),
279 150 : inv_prj_height*(height[2]-(s[2]/s.modulo2())*dotProduct(height,s)));
280 :
281 150 : Value* valuelp=getPntrToComponent("lp");
282 150 : Value* valueld=getPntrToComponent("ld");
283 150 : valuelp->set(dotProduct(s,v)/s.modulo()); // this includes the sign
284 : valueld->set(prj_height);
285 :
286 : // DEBUG
287 : // log.printf(" Dopo: %13.6lf %13.6lf\n",dotProduct(s,v)/s.modulo(),prj_height );
288 :
289 150 : setAtomsDerivatives(valuelp,getNumberOfAtoms()-1,matmul(transpose(Rotation),der_prj));
290 150 : setAtomsDerivatives(valueld,getNumberOfAtoms()-1,matmul(transpose(Rotation),der_height));
291 :
292 150 : Vector der_h;
293 150 : Vector der_l;
294 150 : double weight=1./float(getNumberOfAtoms()-2.);
295 :
296 485850 : for(unsigned iat=0; iat<getNumberOfAtoms()-2; iat++) {
297 : der_h.zero();
298 : der_l.zero();
299 1942800 : for(unsigned a=0; a<3; a++) {
300 5828400 : for(unsigned b=0; b<3; b++) {
301 17485200 : for(unsigned g=0; g<3; g++) {
302 13113900 : der_h[a]+=der_height[b]*drotdpos[b][g][iat][a]*ligand_centered[g];
303 13113900 : der_l[a]+=der_prj[b]*drotdpos[b][g][iat][a]*ligand_centered[g];
304 : }
305 4371300 : der_h[a]-=der_height[b]*Rotation[b][a]*weight;
306 4371300 : der_l[a]-=der_prj[b]*Rotation[b][a]*weight;
307 : }
308 : }
309 485700 : setAtomsDerivatives(valuelp,iat,der_l);
310 485700 : setAtomsDerivatives(valueld,iat,der_h);
311 : }
312 :
313 150 : setBoxDerivativesNoPbc(valuelp);
314 150 : setBoxDerivativesNoPbc(valueld);
315 :
316 300 : }
317 :
318 : }
319 : }
320 :
321 :
322 :
|
