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 FUNNELMOD_COLVAR FUNNEL_PS
26 : /*
27 : FUNNEL_PS implements the Funnel-Metadynamics (FM) technique in PLUMED 2.
28 :
29 : Please read the FM \cite limongelli2013funnel \cite raniolo2020ligand papers 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 \ref FUNNEL 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 \cite kearsley algorithm.
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 reported in
42 : the \ref dists, 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 \ref WHOLEMOLECULES. 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 \ref WHOLEMOLECULES, 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 \ref FUNNEL 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 : \par 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 : \plumedfile
69 : lig: COM ATOMS=2446,2447,2448,2449,2451
70 : fps: FUNNEL_PS REFERENCE=protein.pdb LIGAND=lig POINTS=5.3478,-0.7278,2.4746,7.3785,6.7364,-9.3624
71 : PRINT ARG=fps.lp,fps.ld
72 : \endplumedfile
73 :
74 : It is recommended to add a line to force the reconstruction of the periodic boundary conditions. In the following example,
75 : \ref WHOLEMOLECULES was added to make sure that the ligand was reconstructed together with the protein. The list contains
76 : all the atoms reported in the start.pdb file followed by the ANCHOR atom and the ligand. All atoms should be contained in the
77 : same entity and the correct order.
78 : \plumedfile
79 : 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,
80 : 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,
81 : 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,
82 : 2911,2927,2948,2962,2472,3221,3224,3225,3228,3229,3231,3233,3235,3237
83 : lig: COM ATOMS=3221,3224,3225,3228,3229,3231,3233,3235,3237
84 : fps: FUNNEL_PS LIGAND=lig REFERENCE=start.pdb ANCHOR=2472 POINTS=4.724,5.369,4.069,4.597,5.721,4.343
85 : PRINT ARG=fps.lp,fps.ld
86 : \endplumedfile
87 :
88 : */
89 : //+ENDPLUMEDOC
90 :
91 :
92 : class FUNNEL_PS : public Colvar {
93 :
94 : // Those are arrays that I created for the colvar
95 : std::vector<AtomNumber> ligand_com;
96 : std::vector<AtomNumber> anchor;
97 : std::vector<AtomNumber> numbers;
98 : bool pbc;
99 : PLMD::RMSD alignment;
100 : PLMD::PDB pdb;
101 : bool squared;
102 : private:
103 : vector<double> points;
104 : public:
105 : explicit FUNNEL_PS(const ActionOptions&);
106 : // active methods:
107 : virtual void calculate();
108 : static void registerKeywords(Keywords& keys);
109 : // I need a method in RMSDCoreCalc and these were requested
110 : std::vector<double> align;
111 : std::vector<double> displace;
112 : };
113 :
114 : using namespace std;
115 :
116 : PLUMED_REGISTER_ACTION(FUNNEL_PS,"FUNNEL_PS")
117 :
118 7 : void FUNNEL_PS::registerKeywords(Keywords& keys) {
119 7 : Colvar::registerKeywords( keys );
120 14 : keys.add("compulsory","REFERENCE","a file in pdb format containing the structure you would like to align.");
121 14 : keys.add("atoms","LIGAND","This MUST be a single atom, normally the COM of the ligand");
122 14 : keys.add("atoms","ANCHOR","Closest protein atom to the ligand, picked to avoid pbc problems during the simulation");
123 14 : 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.");
124 14 : keys.addFlag("SQUARED-ROOT",false,"Used to initialize the creation of the alignment variable");
125 14 : keys.addOutputComponent("lp","default","the position along the funnel line");
126 14 : keys.addOutputComponent("ld","default","the distance from the funnel line");
127 7 : }
128 :
129 5 : FUNNEL_PS::FUNNEL_PS(const ActionOptions&ao):
130 : PLUMED_COLVAR_INIT(ao),
131 5 : pbc(true),squared(true)
132 : {
133 : string reference;
134 10 : parse("REFERENCE",reference);
135 : string type;
136 5 : type.assign("OPTIMAL");
137 10 : parseAtomList("LIGAND",ligand_com);
138 5 : if(ligand_com.size()!=1)
139 0 : error("Number of specified atoms should be one, normally the COM of the ligand");
140 5 : parseVector("POINTS",points);
141 5 : bool nopbc=!pbc;
142 5 : parseFlag("NOPBC",nopbc);
143 5 : pbc=!nopbc;
144 : bool sq;
145 5 : parseFlag("SQUARED-ROOT",sq);
146 5 : if (sq) {
147 0 : squared=false;
148 : }
149 5 : parseAtomList("ANCHOR",anchor);
150 5 : checkRead();
151 :
152 : // read everything in ang and transform to nm if we are not in natural units
153 5 : if( !pdb.read(reference,usingNaturalUnits(),0.1/getUnits().getLength()) )
154 0 : error("missing input file " + reference );
155 :
156 : bool remove_com=true;
157 : bool normalize_weights=true;
158 : // here displace is a simple vector of ones
159 5 : align=pdb.getOccupancy();
160 16195 : for(unsigned i=0; i<align.size(); i++) {
161 16190 : align[i]=1.;
162 : } ;
163 5 : displace=pdb.getBeta();
164 16195 : for(unsigned i=0; i<displace.size(); i++) {
165 16190 : displace[i]=1.;
166 : } ;
167 : // reset again to reimpose uniform weights (safe to disable this)
168 5 : alignment.set(align,displace,pdb.getPositions(),type,remove_com,normalize_weights);
169 :
170 :
171 :
172 : // Array with inside both the structure to align and the atom to be aligned
173 5 : numbers=pdb.getAtomNumbers();
174 5 : numbers.push_back(anchor[0]);
175 5 : numbers.push_back(ligand_com[0]);
176 :
177 5 : log.printf(" average from file %s\n",reference.c_str());
178 5 : log.printf(" method for alignment : %s \n",type.c_str() );
179 :
180 5 : if(pbc) log.printf(" using periodic boundary conditions\n");
181 0 : else log.printf(" without periodic boundary conditions\n");
182 :
183 10 : addComponentWithDerivatives("lp");
184 10 : componentIsNotPeriodic("lp");
185 10 : addComponentWithDerivatives("ld");
186 5 : componentIsNotPeriodic("ld");
187 :
188 5 : requestAtoms( numbers );
189 :
190 5 : }
191 :
192 : // calculator
193 150 : void FUNNEL_PS::calculate() {
194 :
195 150 : if(pbc) makeWhole();
196 :
197 150 : Tensor Rotation;
198 : Matrix<std::vector<Vector> > drotdpos(3,3);
199 : std::vector<Vector> alignedpos;
200 : std::vector<Vector> centeredpos;
201 : std::vector<Vector> centeredref;
202 : std::vector<Vector> ddistdpos;
203 :
204 : std::vector<Vector> buffer;
205 :
206 150 : Vector centerreference;
207 150 : Vector centerpositions;
208 :
209 : // Created only to give the correct object to calc_FitElements
210 : std::vector<Vector> sourceAllPositions;
211 : std::vector<Vector> sourcePositions;
212 :
213 : // SourcePositions contains only the coordinates of the COM of the ligand, we need only this point
214 150 : sourceAllPositions=getPositions();
215 150 : sourcePositions=sourceAllPositions;
216 150 : sourcePositions.resize(sourcePositions.size()-2);// just the protein
217 :
218 : // The two points that define the axis : this can be moved in the initialization
219 150 : Vector p1 = VectorGeneric<3>(points[0],points[1],points[2]);
220 150 : Vector p2 = VectorGeneric<3>(points[3],points[4],points[5]);
221 150 : Vector s = p2 - p1;
222 :
223 : // I call the method calc_FitElements that initializes all feature that I need
224 : // except for centerreference that I need to calculate from scratch
225 : // Buffer has no meaning but I had to fulfill the requirements of calc_FitElements
226 150 : double rmsd = alignment.calc_FitElements( sourcePositions, Rotation, drotdpos, buffer, centerpositions, squared);
227 :
228 : // To Plumed developers: it would be interesting to make the functions to calculate centers of mass public or protected
229 150 : centerreference.zero();
230 485850 : for(unsigned i=0; i<pdb.size(); i++) {
231 485700 : centerreference+=pdb.getPositions()[i]*align[i]/align.size();
232 : }
233 :
234 : /*
235 : // I cancelled the additional lines in the library of RMSD.h, thus I am missing the center of the reference
236 : // Creating variable kito to extract only the center of the reference, since no method is calling
237 : // function getReferenceCenter()
238 : PLMD::RMSDCoreData* kito; kito = new RMSDCoreData(align,displace,sourcePositions,pdb->getPositions());
239 : centerreference = kito->getReferenceCenter();
240 : delete kito;
241 : */
242 :
243 : // DEBUG
244 : /* log.printf(" RMSD: %13.6lf\n",rmsd );
245 : log.printf(" cpos: %13.6lf %13.6lf %13.6lf\n",centerpositions[0],centerpositions[1],centerpositions[2] );
246 : log.printf(" Rotation: %13.6lf %13.6lf %13.6lf\n",Rotation[0][0],Rotation[0][1],Rotation[0][2] );
247 : log.printf(" %13.6lf %13.6lf %13.6lf\n",Rotation[1][0],Rotation[1][1],Rotation[1][2] );
248 : log.printf(" %13.6lf %13.6lf %13.6lf\n",Rotation[2][0],Rotation[2][1],Rotation[2][2] );
249 : */
250 :
251 : // the position is Rot(ligand-com_prot_md)+ com_prot_ref
252 150 : Vector ligand_centered =getPositions().back()-centerpositions;
253 150 : Vector ligand_aligned =matmul(Rotation,ligand_centered);
254 150 : Vector v = ligand_aligned +centerreference -p1;
255 :
256 : // DEBUG
257 : // log.printf(" ligando: %13.6lf %13.6lf %13.6lf\n",getPositions().back()[0],getPositions().back()[1],getPositions().back()[2] );
258 :
259 : //Projection vector v onto s
260 :
261 150 : Vector prj = (dotProduct(s,v)/dotProduct(s,s))*s;
262 150 : const double prj_length = prj.modulo() ;
263 : const double inv_prj_length = 1.0/prj_length;
264 :
265 150 : Vector height = v - prj;
266 150 : const double prj_height = height.modulo() ;
267 150 : const double inv_prj_height = 1.0/prj_height;
268 :
269 : // derivative of the prj: only on the com of the ligand
270 150 : Vector der_prj;
271 150 : der_prj=s/s.modulo();
272 :
273 : // derivative of the height: only on the com of the ligand
274 150 : Vector der_height(inv_prj_height*(height[0]-(s[0]/s.modulo2())*dotProduct(height,s)),
275 150 : inv_prj_height*(height[1]-(s[1]/s.modulo2())*dotProduct(height,s)),
276 150 : inv_prj_height*(height[2]-(s[2]/s.modulo2())*dotProduct(height,s)));
277 :
278 150 : Value* valuelp=getPntrToComponent("lp");
279 150 : Value* valueld=getPntrToComponent("ld");
280 150 : valuelp->set(dotProduct(s,v)/s.modulo()); // this includes the sign
281 : valueld->set(prj_height);
282 :
283 : // DEBUG
284 : // log.printf(" Dopo: %13.6lf %13.6lf\n",dotProduct(s,v)/s.modulo(),prj_height );
285 :
286 150 : setAtomsDerivatives(valuelp,getNumberOfAtoms()-1,matmul(transpose(Rotation),der_prj));
287 150 : setAtomsDerivatives(valueld,getNumberOfAtoms()-1,matmul(transpose(Rotation),der_height));
288 :
289 150 : Vector der_h;
290 150 : Vector der_l;
291 150 : double weight=1./float(getNumberOfAtoms()-2.);
292 :
293 485850 : for(unsigned iat=0; iat<getNumberOfAtoms()-2; iat++) {
294 485700 : der_h.zero();
295 485700 : der_l.zero();
296 1942800 : for(unsigned a=0; a<3; a++) {
297 5828400 : for(unsigned b=0; b<3; b++) {
298 17485200 : for(unsigned g=0; g<3; g++) {
299 13113900 : der_h[a]+=der_height[b]*drotdpos[b][g][iat][a]*ligand_centered[g];
300 13113900 : der_l[a]+=der_prj[b]*drotdpos[b][g][iat][a]*ligand_centered[g];
301 : }
302 4371300 : der_h[a]-=der_height[b]*Rotation[b][a]*weight;
303 4371300 : der_l[a]-=der_prj[b]*Rotation[b][a]*weight;
304 : }
305 : }
306 485700 : setAtomsDerivatives(valuelp,iat,der_l);
307 485700 : setAtomsDerivatives(valueld,iat,der_h);
308 : }
309 :
310 150 : setBoxDerivativesNoPbc(valuelp);
311 150 : setBoxDerivativesNoPbc(valueld);
312 :
313 150 : }
314 :
315 : }
316 : }
317 :
318 :
319 :
|
