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 7 : keys.add("compulsory","REFERENCE","a file in pdb format containing the structure you would like to align.");
121 7 : keys.add("atoms","LIGAND","This MUST be a single atom, normally the COM of the ligand");
122 7 : keys.add("atoms","ANCHOR","Closest protein atom to the ligand, picked to avoid pbc problems during the simulation");
123 7 : 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 7 : keys.addFlag("SQUARED-ROOT",false,"Used to initialize the creation of the alignment variable");
125 14 : keys.addOutputComponent("lp","default","scalar","the position along the funnel line");
126 14 : keys.addOutputComponent("ld","default","scalar","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 : string reference;
133 10 : parse("REFERENCE",reference);
134 : string type;
135 5 : type.assign("OPTIMAL");
136 10 : parseAtomList("LIGAND",ligand_com);
137 5 : if(ligand_com.size()!=1) {
138 0 : error("Number of specified atoms should be one, normally the COM of the ligand");
139 : }
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 :
157 : bool remove_com=true;
158 : bool normalize_weights=true;
159 : // here displace is a simple vector of ones
160 5 : align=pdb.getOccupancy();
161 16195 : for(unsigned i=0; i<align.size(); i++) {
162 16190 : align[i]=1.;
163 : } ;
164 5 : displace=pdb.getBeta();
165 16195 : for(unsigned i=0; i<displace.size(); i++) {
166 16190 : displace[i]=1.;
167 : } ;
168 : // reset again to reimpose uniform weights (safe to disable this)
169 5 : alignment.set(align,displace,pdb.getPositions(),type,remove_com,normalize_weights);
170 :
171 :
172 :
173 : // Array with inside both the structure to align and the atom to be aligned
174 5 : numbers=pdb.getAtomNumbers();
175 5 : numbers.push_back(anchor[0]);
176 5 : numbers.push_back(ligand_com[0]);
177 :
178 5 : log.printf(" average from file %s\n",reference.c_str());
179 5 : log.printf(" method for alignment : %s \n",type.c_str() );
180 :
181 5 : if(pbc) {
182 5 : log.printf(" using periodic boundary conditions\n");
183 : } else {
184 0 : log.printf(" without periodic boundary conditions\n");
185 : }
186 :
187 10 : addComponentWithDerivatives("lp");
188 10 : componentIsNotPeriodic("lp");
189 10 : addComponentWithDerivatives("ld");
190 5 : componentIsNotPeriodic("ld");
191 :
192 5 : requestAtoms( numbers );
193 :
194 5 : }
195 :
196 : // calculator
197 150 : void FUNNEL_PS::calculate() {
198 :
199 150 : if(pbc) {
200 150 : makeWhole();
201 : }
202 :
203 150 : Tensor Rotation;
204 : Matrix<std::vector<Vector> > drotdpos(3,3);
205 : std::vector<Vector> alignedpos;
206 : std::vector<Vector> centeredpos;
207 : std::vector<Vector> centeredref;
208 : std::vector<Vector> ddistdpos;
209 :
210 : std::vector<Vector> buffer;
211 :
212 150 : Vector centerreference;
213 150 : Vector centerpositions;
214 :
215 : // Created only to give the correct object to calc_FitElements
216 : std::vector<Vector> sourceAllPositions;
217 : std::vector<Vector> sourcePositions;
218 :
219 : // SourcePositions contains only the coordinates of the COM of the ligand, we need only this point
220 150 : sourceAllPositions=getPositions();
221 150 : sourcePositions=sourceAllPositions;
222 150 : sourcePositions.resize(sourcePositions.size()-2);// just the protein
223 :
224 : // The two points that define the axis : this can be moved in the initialization
225 150 : Vector p1 = VectorGeneric<3>(points[0],points[1],points[2]);
226 150 : Vector p2 = VectorGeneric<3>(points[3],points[4],points[5]);
227 150 : Vector s = p2 - p1;
228 :
229 : // I call the method calc_FitElements that initializes all feature that I need
230 : // except for centerreference that I need to calculate from scratch
231 : // Buffer has no meaning but I had to fulfill the requirements of calc_FitElements
232 150 : double rmsd = alignment.calc_FitElements( sourcePositions, Rotation, drotdpos, buffer, centerpositions, squared);
233 :
234 : // To Plumed developers: it would be interesting to make the functions to calculate centers of mass public or protected
235 150 : centerreference.zero();
236 485850 : for(unsigned i=0; i<pdb.size(); i++) {
237 485700 : centerreference+=pdb.getPositions()[i]*align[i]/align.size();
238 : }
239 :
240 : /*
241 : // I cancelled the additional lines in the library of RMSD.h, thus I am missing the center of the reference
242 : // Creating variable kito to extract only the center of the reference, since no method is calling
243 : // function getReferenceCenter()
244 : PLMD::RMSDCoreData* kito; kito = new RMSDCoreData(align,displace,sourcePositions,pdb->getPositions());
245 : centerreference = kito->getReferenceCenter();
246 : delete kito;
247 : */
248 :
249 : // DEBUG
250 : /* log.printf(" RMSD: %13.6lf\n",rmsd );
251 : log.printf(" cpos: %13.6lf %13.6lf %13.6lf\n",centerpositions[0],centerpositions[1],centerpositions[2] );
252 : log.printf(" Rotation: %13.6lf %13.6lf %13.6lf\n",Rotation[0][0],Rotation[0][1],Rotation[0][2] );
253 : log.printf(" %13.6lf %13.6lf %13.6lf\n",Rotation[1][0],Rotation[1][1],Rotation[1][2] );
254 : log.printf(" %13.6lf %13.6lf %13.6lf\n",Rotation[2][0],Rotation[2][1],Rotation[2][2] );
255 : */
256 :
257 : // the position is Rot(ligand-com_prot_md)+ com_prot_ref
258 150 : Vector ligand_centered =getPositions().back()-centerpositions;
259 150 : Vector ligand_aligned =matmul(Rotation,ligand_centered);
260 150 : Vector v = ligand_aligned +centerreference -p1;
261 :
262 : // DEBUG
263 : // log.printf(" ligando: %13.6lf %13.6lf %13.6lf\n",getPositions().back()[0],getPositions().back()[1],getPositions().back()[2] );
264 :
265 : //Projection vector v onto s
266 :
267 150 : Vector prj = (dotProduct(s,v)/dotProduct(s,s))*s;
268 150 : const double prj_length = prj.modulo() ;
269 : const double inv_prj_length = 1.0/prj_length;
270 :
271 150 : Vector height = v - prj;
272 150 : const double prj_height = height.modulo() ;
273 150 : const double inv_prj_height = 1.0/prj_height;
274 :
275 : // derivative of the prj: only on the com of the ligand
276 150 : Vector der_prj;
277 150 : der_prj=s/s.modulo();
278 :
279 : // derivative of the height: only on the com of the ligand
280 150 : Vector der_height(inv_prj_height*(height[0]-(s[0]/s.modulo2())*dotProduct(height,s)),
281 150 : inv_prj_height*(height[1]-(s[1]/s.modulo2())*dotProduct(height,s)),
282 150 : inv_prj_height*(height[2]-(s[2]/s.modulo2())*dotProduct(height,s)));
283 :
284 150 : Value* valuelp=getPntrToComponent("lp");
285 150 : Value* valueld=getPntrToComponent("ld");
286 150 : valuelp->set(dotProduct(s,v)/s.modulo()); // this includes the sign
287 : valueld->set(prj_height);
288 :
289 : // DEBUG
290 : // log.printf(" Dopo: %13.6lf %13.6lf\n",dotProduct(s,v)/s.modulo(),prj_height );
291 :
292 150 : setAtomsDerivatives(valuelp,getNumberOfAtoms()-1,matmul(transpose(Rotation),der_prj));
293 150 : setAtomsDerivatives(valueld,getNumberOfAtoms()-1,matmul(transpose(Rotation),der_height));
294 :
295 150 : Vector der_h;
296 150 : Vector der_l;
297 150 : double weight=1./float(getNumberOfAtoms()-2.);
298 :
299 485850 : for(unsigned iat=0; iat<getNumberOfAtoms()-2; iat++) {
300 485700 : der_h.zero();
301 485700 : der_l.zero();
302 1942800 : for(unsigned a=0; a<3; a++) {
303 5828400 : for(unsigned b=0; b<3; b++) {
304 17485200 : for(unsigned g=0; g<3; g++) {
305 13113900 : der_h[a]+=der_height[b]*drotdpos[b][g][iat][a]*ligand_centered[g];
306 13113900 : der_l[a]+=der_prj[b]*drotdpos[b][g][iat][a]*ligand_centered[g];
307 : }
308 4371300 : der_h[a]-=der_height[b]*Rotation[b][a]*weight;
309 4371300 : der_l[a]-=der_prj[b]*Rotation[b][a]*weight;
310 : }
311 : }
312 485700 : setAtomsDerivatives(valuelp,iat,der_l);
313 485700 : setAtomsDerivatives(valueld,iat,der_h);
314 : }
315 :
316 150 : setBoxDerivativesNoPbc(valuelp);
317 150 : setBoxDerivativesNoPbc(valueld);
318 :
319 150 : }
320 :
321 : }
322 : }
323 :
324 :
325 :
|
