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