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 "colvar/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 "core/PlumedMain.h"
18 : #include "core/Atoms.h"
19 : #include <iostream>
20 : #include "tools/RMSD.h"
21 :
22 : using namespace std;
23 :
24 : namespace PLMD {
25 : namespace funnel {
26 :
27 : //+PLUMEDOC FUNNELMOD_COLVAR FUNNEL_PS
28 : /*
29 : FUNNEL_PS implements the Funnel-Metadynamics (FM) technique in PLUMED 2.
30 :
31 : Please read the FM \cite limongelli2013funnel \cite raniolo2020ligand papers to better understand the notions hereby reported.
32 :
33 : This colvar evaluates the position of a ligand of interest with respect to a given line, built from the two points
34 : A and B, and should be used together with the \ref FUNNEL bias.
35 : The constructed line represents the axis of the funnel-shape restraint potential, which should be placed so
36 : as to include the portion of a macromolecule (i.e., protein, DNA, etc.) that should be explored.
37 : Since it is important that the position of the line is updated based on the motion of the macromolecule during
38 : the simulation, this colvar incorporates an alignment method. The latter uses the TYPE=OPTIMAL option to remove
39 : motions due to rotation and translation of the macromolecule with respect to a reference structure, which is
40 : provided by the user. In order to accomplish the task, an optimal alignment matrix is calculated using the
41 : Kearsley \cite kearsley algorithm.
42 : The reference structure should be given as a pdb file, containing only the atoms of the macromolecule or a
43 : selection of them (e.g., the protein CA atoms of secondary structures). In contrast to the methods reported in
44 : the \ref dists, the values reported in the occupancy and beta-factor columns of the pdb file are not important
45 : since they will be overwritten and replaced by the value 1.00 during the procedure. It is important to understand
46 : that all atoms in the file will be used for the alignment, even if they display 0.00 in the occupancy column.
47 :
48 : The ligand can be represented by one single atom or the center of mass (COM) of a group of atoms that should be
49 : provided by the user.
50 :
51 : By default FUNNEL_PS is computed taking into account periodic boundary conditions. Since PLUMED 2.5, molecules are
52 : rebuilt using a procedure that is equivalent to that done in \ref WHOLEMOLECULES. We note that this action is local
53 : to this colvar, thus it does not modify the coordinates stored in PLUMED. Moreover, FUNNEL_PS requires an ANCHOR atom
54 : to be specified in order to facilitate the reconstruction of periodic boundary conditions. This atom should be the
55 : closest macromolecule's atom to the ligand and it should reduce the risk of ligand "warping" in the simulation box.
56 : Nevertheless, we highly recommend to add to the PLUMED input file a custom line of \ref WHOLEMOLECULES, in order to
57 : be sure of reconstructing the ligand together with the macromolecule (please look the examples). In this case, the user
58 : can use the NOPBC flag to turn off the internal periodic boundary condition reconstruction.
59 :
60 : FUNNEL_PS is divided in two components (fps.lp and fps.ld) which evaluate the projection of the ligand along the funnel line
61 : and the distance from it, respectively. The values attributed to these two components are then used together with the
62 : potential file created by the \ref FUNNEL bias to define if the ligand is within or not in the funnel-shape restraint
63 : potential. In the latter case, the potential will force the ligand to enter within the funnel boundaries.
64 :
65 : \par Examples
66 :
67 : The following input tells plumed to print the FUNNEL_PS components for the COM of a ligand. The inputs are a reference structure,
68 : which is the structure used for the alignment, the COM of the molecule we want to track, and 2 points in the Cartesian space
69 : (i.e., x, y, and z) to draw the line representing the funnel axis.
70 : \plumedfile
71 : lig: COM ATOMS=2446,2447,2448,2449,2451
72 : fps: FUNNEL_PS REFERENCE=protein.pdb LIGAND=lig POINTS=5.3478,-0.7278,2.4746,7.3785,6.7364,-9.3624
73 : PRINT ARG=fps.lp,fps.ld
74 : \endplumedfile
75 :
76 : It is recommended to add a line to force the reconstruction of the periodic boundary conditions. In the following example,
77 : \ref WHOLEMOLECULES was added to make sure that the ligand was reconstructed together with the protein. The list contains
78 : all the atoms reported in the start.pdb file followed by the ANCHOR atom and the ligand. All atoms should be contained in the
79 : same entity and the correct order.
80 : \plumedfile
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 : \endplumedfile
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 : // It is written no more desctructors, but an expert said it's necessary for imported variables (pdb and alignment) or else memory leak
115 : ~FUNNEL_PS();
116 : };
117 :
118 : using namespace std;
119 :
120 10429 : PLUMED_REGISTER_ACTION(FUNNEL_PS,"FUNNEL_PS")
121 :
122 6 : void FUNNEL_PS::registerKeywords(Keywords& keys) {
123 6 : Colvar::registerKeywords( keys );
124 12 : keys.add("compulsory","REFERENCE","a file in pdb format containing the structure you would like to align.");
125 12 : keys.add("atoms","LIGAND","This MUST be a single atom, normally the COM of the ligand");
126 12 : keys.add("atoms","ANCHOR","Closest protein atom to the ligand, picked to avoid pbc problems during the simulation");
127 12 : 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.");
128 12 : keys.addFlag("SQUARED-ROOT",false,"Used to initialize the creation of the alignment variable");
129 12 : keys.addOutputComponent("lp","default","the position along the funnel line");
130 12 : keys.addOutputComponent("ld","default","the distance from the funnel line");
131 6 : }
132 :
133 5 : FUNNEL_PS::FUNNEL_PS(const ActionOptions&ao):
134 : PLUMED_COLVAR_INIT(ao),
135 5 : pbc(true),squared(true)
136 : {
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 5 : parseVector("POINTS",points);
145 5 : bool nopbc=!pbc;
146 5 : parseFlag("NOPBC",nopbc);
147 5 : pbc=!nopbc;
148 : bool sq;
149 5 : parseFlag("SQUARED-ROOT",sq);
150 5 : if (sq) {
151 0 : squared=false;
152 : }
153 5 : parseAtomList("ANCHOR",anchor);
154 5 : checkRead();
155 :
156 5 : pdb = new PDB();
157 :
158 : // read everything in ang and transform to nm if we are not in natural units
159 10 : if( !pdb->read(reference,plumed.getAtoms().usingNaturalUnits(),0.1/ActionAtomistic::atoms.getUnits().getLength()) )
160 0 : error("missing input file " + reference );
161 :
162 5 : alignment = new RMSD();
163 :
164 : bool remove_com=true;
165 : bool normalize_weights=true;
166 : // here displace is a simple vector of ones
167 5 : align=pdb->getOccupancy();
168 16195 : for(unsigned i=0; i<align.size(); i++) {
169 16190 : align[i]=1.;
170 : } ;
171 5 : displace=pdb->getBeta();
172 16195 : for(unsigned i=0; i<displace.size(); i++) {
173 16190 : displace[i]=1.;
174 : } ;
175 : // reset again to reimpose uniform weights (safe to disable this)
176 5 : alignment->set(align,displace,pdb->getPositions(),type,remove_com,normalize_weights);
177 :
178 :
179 :
180 : // Array with inside both the structure to align and the atom to be aligned
181 5 : numbers=pdb->getAtomNumbers();
182 5 : numbers.push_back(anchor[0]);
183 5 : numbers.push_back(ligand_com[0]);
184 :
185 5 : log.printf(" average from file %s\n",reference.c_str());
186 5 : log.printf(" method for alignment : %s \n",type.c_str() );
187 :
188 5 : if(pbc) log.printf(" using periodic boundary conditions\n");
189 0 : else log.printf(" without periodic boundary conditions\n");
190 :
191 5 : addComponentWithDerivatives("lp");
192 5 : componentIsNotPeriodic("lp");
193 5 : addComponentWithDerivatives("ld");
194 5 : componentIsNotPeriodic("ld");
195 :
196 5 : requestAtoms( numbers );
197 :
198 5 : }
199 :
200 10 : FUNNEL_PS::~FUNNEL_PS() {
201 5 : delete alignment;
202 5 : delete pdb;
203 10 : }
204 :
205 : // calculator
206 150 : void FUNNEL_PS::calculate() {
207 :
208 150 : if(pbc) makeWhole();
209 :
210 150 : Tensor Rotation;
211 : Matrix<std::vector<Vector> > drotdpos(3,3);
212 : std::vector<Vector> alignedpos;
213 : std::vector<Vector> centeredpos;
214 : std::vector<Vector> centeredref;
215 : std::vector<Vector> ddistdpos;
216 :
217 : std::vector<Vector> buffer;
218 :
219 150 : Vector centerreference;
220 150 : Vector centerpositions;
221 :
222 : // Created only to give the correct object to calc_FitElements
223 : std::vector<Vector> sourceAllPositions;
224 : std::vector<Vector> sourcePositions;
225 :
226 : // SourcePositions contains only the coordinates of the COM of the ligand, we need only this point
227 150 : sourceAllPositions=getPositions();
228 150 : sourcePositions=sourceAllPositions;
229 150 : sourcePositions.resize(sourcePositions.size()-2);// just the protein
230 :
231 : // The two points that define the axis : this can be moved in the initialization
232 150 : Vector p1 = VectorGeneric<3>(points[0],points[1],points[2]);
233 150 : Vector p2 = VectorGeneric<3>(points[3],points[4],points[5]);
234 150 : Vector s = p2 - p1;
235 :
236 : // I call the method calc_FitElements that initializes all feature that I need
237 : // except for centerreference that I need to calculate from scratch
238 : // Buffer has no meaning but I had to fulfill the requirements of calc_FitElements
239 150 : double rmsd = alignment->calc_FitElements( sourcePositions, Rotation, drotdpos, buffer, centerpositions, squared);
240 :
241 : // To Plumed developers: it would be interesting to make the functions to calculate centers of mass public or protected
242 150 : centerreference.zero();
243 485850 : for(unsigned i=0; i<pdb->size(); i++) {
244 485700 : centerreference+=pdb->getPositions()[i]*align[i]/align.size();
245 : }
246 :
247 : /*
248 : // I cancelled the additional lines in the library of RMSD.h, thus I am missing the center of the reference
249 : // Creating variable kito to extract only the center of the reference, since no method is calling
250 : // function getReferenceCenter()
251 : PLMD::RMSDCoreData* kito; kito = new RMSDCoreData(align,displace,sourcePositions,pdb->getPositions());
252 : centerreference = kito->getReferenceCenter();
253 : delete kito;
254 : */
255 :
256 : // DEBUG
257 : /* log.printf(" RMSD: %13.6lf\n",rmsd );
258 : log.printf(" cpos: %13.6lf %13.6lf %13.6lf\n",centerpositions[0],centerpositions[1],centerpositions[2] );
259 : log.printf(" Rotation: %13.6lf %13.6lf %13.6lf\n",Rotation[0][0],Rotation[0][1],Rotation[0][2] );
260 : log.printf(" %13.6lf %13.6lf %13.6lf\n",Rotation[1][0],Rotation[1][1],Rotation[1][2] );
261 : log.printf(" %13.6lf %13.6lf %13.6lf\n",Rotation[2][0],Rotation[2][1],Rotation[2][2] );
262 : */
263 :
264 : // the position is Rot(ligand-com_prot_md)+ com_prot_ref
265 150 : Vector ligand_centered =getPositions().back()-centerpositions;
266 150 : Vector ligand_aligned =matmul(Rotation,ligand_centered);
267 150 : Vector v = ligand_aligned +centerreference -p1;
268 :
269 : // DEBUG
270 : // log.printf(" ligando: %13.6lf %13.6lf %13.6lf\n",getPositions().back()[0],getPositions().back()[1],getPositions().back()[2] );
271 :
272 : //Projection vector v onto s
273 :
274 150 : Vector prj = (dotProduct(s,v)/dotProduct(s,s))*s;
275 150 : const double prj_length = prj.modulo() ;
276 : const double inv_prj_length = 1.0/prj_length;
277 :
278 150 : Vector height = v - prj;
279 150 : const double prj_height = height.modulo() ;
280 150 : const double inv_prj_height = 1.0/prj_height;
281 :
282 : // derivative of the prj: only on the com of the ligand
283 150 : Vector der_prj;
284 150 : der_prj=s/s.modulo();
285 :
286 : // derivative of the height: only on the com of the ligand
287 150 : Vector der_height(inv_prj_height*(height[0]-(s[0]/s.modulo2())*dotProduct(height,s)),
288 150 : inv_prj_height*(height[1]-(s[1]/s.modulo2())*dotProduct(height,s)),
289 150 : inv_prj_height*(height[2]-(s[2]/s.modulo2())*dotProduct(height,s)));
290 :
291 150 : Value* valuelp=getPntrToComponent("lp");
292 150 : Value* valueld=getPntrToComponent("ld");
293 150 : valuelp->set(dotProduct(s,v)/s.modulo()); // this includes the sign
294 : valueld->set(prj_height);
295 :
296 : // DEBUG
297 : // log.printf(" Dopo: %13.6lf %13.6lf\n",dotProduct(s,v)/s.modulo(),prj_height );
298 :
299 150 : setAtomsDerivatives(valuelp,getNumberOfAtoms()-1,matmul(transpose(Rotation),der_prj));
300 150 : setAtomsDerivatives(valueld,getNumberOfAtoms()-1,matmul(transpose(Rotation),der_height));
301 :
302 150 : Vector der_h;
303 150 : Vector der_l;
304 150 : double weight=1./float(getNumberOfAtoms()-2.);
305 :
306 485850 : for(unsigned iat=0; iat<getNumberOfAtoms()-2; iat++) {
307 485700 : der_h.zero();
308 485700 : der_l.zero();
309 1942800 : for(unsigned a=0; a<3; a++) {
310 5828400 : for(unsigned b=0; b<3; b++) {
311 17485200 : for(unsigned g=0; g<3; g++) {
312 13113900 : der_h[a]+=der_height[b]*drotdpos[b][g][iat][a]*ligand_centered[g];
313 13113900 : der_l[a]+=der_prj[b]*drotdpos[b][g][iat][a]*ligand_centered[g];
314 : }
315 4371300 : der_h[a]-=der_height[b]*Rotation[b][a]*weight;
316 4371300 : der_l[a]-=der_prj[b]*Rotation[b][a]*weight;
317 : }
318 : }
319 485700 : setAtomsDerivatives(valuelp,iat,der_l);
320 485700 : setAtomsDerivatives(valueld,iat,der_h);
321 : }
322 :
323 150 : setBoxDerivativesNoPbc(valuelp);
324 150 : setBoxDerivativesNoPbc(valueld);
325 :
326 150 : }
327 :
328 : }
329 : }
330 :
331 :
332 :
|
