Line data Source code
1 : /* +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
2 : Copyright (c) 2024 by Glen Hocky, New York University on behalf of authors
3 :
4 : The sizeshape module is free software: you can redistribute it and/or modify
5 : it under the terms of the GNU Lesser General Public License as published by
6 : the Free Software Foundation, either version 3 of the License, or
7 : (at your option) any later version.
8 :
9 : The sizeshape module is distributed in the hope that it will be useful,
10 : but WITHOUT ANY WARRANTY; without even the implied warranty of
11 : MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
12 : GNU Lesser General Public License for more details.
13 :
14 : You should have received a copy of the GNU Lesser General Public License
15 : along with plumed. If not, see <http://www.gnu.org/licenses/>.
16 : +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ */
17 : #include "colvar/Colvar.h"
18 : #include "core/ActionRegister.h"
19 : #include "tools/Pbc.h"
20 : #include "tools/File.h" // Input and output from files
21 : #include "tools/Matrix.h" // Linear Algebra operations
22 : #include <sstream>
23 : #include <cmath>
24 : #include "tools/Communicator.h" // All the MPI related stuffs
25 :
26 : namespace PLMD {
27 : namespace sizeshape {
28 :
29 : //+PLUMEDOC sizeshapeMOD_COLVAR SIZESHAPE_POSITION_LINEAR_PROJ
30 : /*
31 : Calculates a linear projection in the space of a given reference configurational distribution in size-and-shape space.
32 :
33 : This method is described in \cite Sasmal-poslda-2023.
34 :
35 : The linear projection is given by:
36 : \f[
37 : l(\mathbf{x}) = \mathbf{v}\cdot(\mathbf{R}\cdot(\mathbf{x}(t) - \vec{{\zeta}}(t)) - \mathbf{\mu}),
38 : \f]
39 : Where \f$\mathbf{v}\f$ is a vector of linear coefficients, \f$\mathbf{x}(t)\f$ is the configuration at time t, \f$\vec{\zeta}(t)\f$ is the difference in the geometric mean of the current configuration and that of the reference configuration \f$\mathbf{\mu}\f$. \f$\vec{\zeta}(t) = \frac{1}{N} \sum_{i=1}^N \vec{x_{i}}(t) - \frac{1}{N} \sum_{i=1}^N \vec{\mu_{i}}(t)\f$, for N atoms.
40 :
41 : \f$\mathbf{R}\f$ is an optimal rotation matrix that minimizes the Mahalanobis distance between the current configuration and reference. \f$\mathbf{R}\f$ is obtained by using Kabsch algorithm within the code. The Mahalanobis distance is given as:
42 :
43 : \f[
44 : d(\mathbf{x}, \mathbf{\mu}, \mathbf{\Sigma}) = \sqrt{(\mathbf{x}-\mathbf{\mu})^T \mathbf{\Sigma}^{-1} (\mathbf{x}-\mathbf{\mu})}
45 : \f]
46 :
47 : where, \f$\mathbf{\Sigma}^{-1}\f$ is the \f$N\times N\f$ precision matrix. See also \ref POSITION_MAHA_DIST for information about calculating Mahalanobis distance in size-and-shape space.
48 :
49 : Size-and-shape Gaussian Mixture Model (shapeGMM) \cite Heidi-shapeGMM-2022 is a probabilistic clustering technique that is used to perform structural clusteing on ensemble of molecular configurations and to obtain reference \f$(\mathbf{\mu})\f$ and precision \f$(\mathbf{\Sigma}^{-1})\f$ corresponding to each of the cluster centers. Please chcek out <a href="https://github.com/mccullaghlab/shapeGMMTorch">shapeGMMTorch-GitHub</a> and <a href="https://pypi.org/project/shapeGMMTorch/"> shapeGMMTorch-PyPI</a> for examples and informations on preforming shapeGMM clustering.
50 :
51 : \par Examples
52 : In the following example, a group is defined with atom indices of selected atoms and then linear projection is calculated for the given reference, precision and coefficients. Each file is a space separated list of 3N floating point numbers.
53 :
54 : \plumedfile
55 : UNITS LENGTH=A TIME=ps ENERGY=kcal/mol
56 : GROUP ATOMS=18,20,22,31,33,35,44,46,48,57,59,61,70,72,74,83,85,87,96,98,100,109,111 LABEL=ga_list
57 : #SETTINGS AUXFILE=regtest/sizeshape/rt-sizeshape/global_avg.txt
58 : #SETTINGS AUXFILE=regtest/sizeshape/rt-sizeshape/global_precision.txt
59 : #SETTINGS AUXFILE=regtest/sizeshape/rt-sizeshape/ld1_scalings.txt
60 : proj: SIZESHAPE_POSITION_LINEAR_PROJ REFERENCE=global_avg.txt PRECISION=global_precision.txt COEFFS=ld1_scalings.txt GROUP=ga_list
61 : PRINT ARG=proj STRIDE=1 FILE=COLVAR FMT=%8.8f
62 : \endplumedfile
63 :
64 : */
65 : //+ENDPLUMEDOC
66 :
67 :
68 : class position_linear_proj : public Colvar {
69 :
70 : private:
71 : bool pbc, serial;
72 : std::string prec_f_name; // precision file name
73 : std::string ref_f_name; // reference file name
74 : std::string coeffs_f_name; // file containing linear coeffs
75 : IFile in_; // create an object of class IFile
76 : Log out_;
77 : Matrix<double> ref_str; // coords of reference
78 : Matrix<double> mobile_str; // coords of mobile
79 : Matrix<double> prec; // precision data
80 : Matrix<double> rotation;
81 : std::vector<double> linear_coeffs; // Linear Coefficients
82 : Matrix<double> derv_numeric;
83 : void readinputs(); // reads the input data
84 : double proj; // projection value
85 : std::vector<AtomNumber> atom_list; // list of atoms
86 : const double SMALL = 1.0E-30;
87 : const double delta = 0.00001;
88 : public:
89 : static void registerKeywords( Keywords& keys );
90 : explicit position_linear_proj(const ActionOptions&);
91 : double determinant(int n, const std::vector<std::vector<double>>* B);
92 : void kabsch_rot_mat(); // gives rotation matrix
93 : double cal_position_linear_proj(); // calculates the linear projection
94 : void numeric_grad(); // calculates the numeric gradient
95 : // active methods:
96 : void calculate() override;
97 : };
98 :
99 : PLUMED_REGISTER_ACTION(position_linear_proj, "SIZESHAPE_POSITION_LINEAR_PROJ")
100 :
101 : // register keywords function
102 7 : void position_linear_proj::registerKeywords( Keywords& keys ) {
103 7 : Colvar::registerKeywords( keys );
104 7 : keys.add("compulsory", "PRECISION", "Precision Matrix (inverse of covariance)." );
105 7 : keys.add("compulsory", "REFERENCE", "Coordinates of the reference structure.");
106 7 : keys.add("atoms","GROUP","Group of atoms being used");
107 7 : keys.add("compulsory", "COEFFS", "Vector of linear coefficients.");
108 7 : keys.addFlag("SERIAL",false,"Perform the calculation in serial, for debug purposes only.");
109 14 : keys.setValueDescription("scalar","the linear projection");
110 7 : }
111 :
112 : // constructor function
113 5 : position_linear_proj::position_linear_proj(const ActionOptions&ao):
114 : PLUMED_COLVAR_INIT(ao),
115 5 : pbc(true),
116 5 : serial(false),
117 5 : proj(0),
118 10 : prec_f_name(""),
119 5 : ref_f_name(""),
120 15 : coeffs_f_name("") { // Note! no comma here in the last line.
121 5 : parseFlag("SERIAL",serial);
122 5 : parseAtomList("GROUP",atom_list);
123 5 : parse("REFERENCE", ref_f_name);
124 5 : parse("PRECISION", prec_f_name);
125 5 : parse("COEFFS", coeffs_f_name);
126 5 : bool nopbc=!pbc;
127 5 : parseFlag("NOPBC",nopbc);
128 5 : pbc=!nopbc;
129 :
130 5 : checkRead();
131 :
132 5 : log.printf(" of %u atoms\n",static_cast<unsigned>(atom_list.size()));
133 120 : for(unsigned int i=0; i<atom_list.size(); ++i) {
134 115 : log.printf(" %d", atom_list[i].serial());
135 : }
136 :
137 5 : if(pbc) {
138 5 : log.printf("\n using periodic boundary conditions\n");
139 : } else {
140 0 : log.printf("\n without periodic boundary conditions\n");
141 : }
142 :
143 5 : addValueWithDerivatives();
144 5 : setNotPeriodic();
145 :
146 5 : requestAtoms(atom_list);
147 :
148 : // call the readinputs() function here
149 5 : readinputs();
150 :
151 5 : }
152 :
153 : // read inputs function
154 5 : void position_linear_proj::readinputs() {
155 : unsigned N=getNumberOfAtoms();
156 : // read ref coords
157 5 : in_.open(ref_f_name);
158 :
159 : ref_str.resize(N,3);
160 : prec.resize(N,N);
161 : derv_numeric.resize(N,3);
162 :
163 : std::string line_, val_;
164 : unsigned c_=0;
165 :
166 120 : while (c_ < N) {
167 115 : in_.getline(line_);
168 : std::vector<std::string> items_;
169 115 : std::stringstream check_(line_);
170 :
171 460 : while(std::getline(check_, val_, ' ')) {
172 345 : items_.push_back(val_);
173 : }
174 460 : for(int i=0; i<3; ++i) {
175 345 : ref_str(c_,i) = std::stold(items_[i]);
176 : }
177 115 : c_ += 1;
178 115 : }
179 5 : in_.close();
180 :
181 : //read precision
182 5 : in_.open(prec_f_name);
183 :
184 : std::string line, val;
185 : unsigned int c = 0;
186 :
187 120 : while(c < N) {
188 115 : in_.getline(line);
189 :
190 : // vector for storing the objects
191 : std::vector<std::string> items;
192 :
193 : // stringstream helps to treat a string like an ifstream!
194 115 : std::stringstream check(line);
195 :
196 2760 : while (std::getline(check, val, ' ')) {
197 2645 : items.push_back(val);
198 : }
199 :
200 2760 : for(unsigned int i=0; i<N; ++i) {
201 2645 : prec(c, i) = std::stold(items[i]);
202 : }
203 :
204 115 : c += 1;
205 :
206 115 : }
207 5 : in_.close();
208 :
209 : // read in the linear coeffs
210 5 : in_.open(coeffs_f_name);
211 : unsigned n_=0;
212 : std::string l_;
213 350 : while (n_ < N*3) {
214 345 : in_.getline(l_);
215 345 : linear_coeffs.push_back(std::stod(l_));
216 345 : n_ += 1;
217 : }
218 5 : linear_coeffs.resize(N*3);
219 :
220 5 : in_.close();
221 :
222 5 : }
223 :
224 :
225 :
226 1430 : double position_linear_proj::determinant(int n, const std::vector<std::vector<double>>* B) {
227 :
228 1430 : std::vector<std::vector<double>> A(n, std::vector<double>(n, 0));
229 : // make a copy first!
230 5720 : for(int i=0; i<n; ++i) {
231 17160 : for(int j=0; j<n; ++j) {
232 12870 : A[i][j] = (*B)[i][j];
233 : }
234 : }
235 :
236 :
237 : // It calculates determinant of a matrix using partial pivoting.
238 :
239 : double det = 1;
240 :
241 : // Row operations for i = 0, ,,,, n - 2 (n-1 not needed)
242 4290 : for ( int i = 0; i < n - 1; i++ ) {
243 : // Partial pivot: find row r below with largest element in column i
244 : int r = i;
245 2860 : double maxA = std::abs( A[i][i] );
246 7150 : for ( int k = i + 1; k < n; k++ ) {
247 4290 : double val = std::abs( A[k][i] );
248 4290 : if ( val > maxA ) {
249 : r = k;
250 : maxA = val;
251 : }
252 : }
253 2860 : if ( r != i ) {
254 10010 : for ( int j = i; j < n; j++ ) {
255 7150 : std::swap( A[i][j], A[r][j] );
256 : }
257 2860 : det = -det;
258 : }
259 :
260 : // Row operations to make upper-triangular
261 2860 : double pivot = A[i][i];
262 2860 : if (std::abs( pivot ) < SMALL ) {
263 : return 0.0; // Singular matrix
264 : }
265 :
266 7150 : for ( int r = i + 1; r < n; r++ ) { // On lower rows
267 4290 : double multiple = A[r][i] / pivot; // Multiple of row i to clear element in ith column
268 15730 : for ( int j = i; j < n; j++ ) {
269 11440 : A[r][j] -= multiple * A[i][j];
270 : }
271 : }
272 2860 : det *= pivot; // Determinant is product of diagonal
273 : }
274 :
275 1430 : det *= A[n-1][n-1];
276 :
277 1430 : return det;
278 1430 : }
279 :
280 : // kabsch rotation
281 715 : void position_linear_proj::kabsch_rot_mat() {
282 :
283 : unsigned N=getNumberOfAtoms();
284 :
285 : Matrix<double> mobile_str_T(3,N);
286 : Matrix<double> prec_dot_ref_str(N,3);
287 : Matrix<double> correlation(3,3);
288 :
289 :
290 715 : transpose(mobile_str, mobile_str_T);
291 715 : mult(prec, ref_str, prec_dot_ref_str);
292 715 : mult(mobile_str_T, prec_dot_ref_str, correlation);
293 :
294 :
295 715 : int rw = correlation.nrows();
296 715 : int cl = correlation.ncols();
297 715 : int sz = rw*cl;
298 :
299 : // SVD part (taking from plu2/src/tools/Matrix.h: pseudoInvert function)
300 :
301 715 : std::vector<double> da(sz);
302 : unsigned k=0;
303 :
304 : // Transfer the matrix to the local array
305 2860 : for (int i=0; i<cl; ++i)
306 8580 : for (int j=0; j<rw; ++j) {
307 6435 : da[k++]=static_cast<double>( correlation(j,i) ); // note! its [j][i] not [i][j]
308 : }
309 :
310 715 : int nsv, info, nrows=rw, ncols=cl;
311 : if(rw>cl) {
312 : nsv=cl;
313 : } else {
314 : nsv=rw;
315 : }
316 :
317 : // Create some containers for stuff from single value decomposition
318 715 : std::vector<double> S(nsv);
319 715 : std::vector<double> U(nrows*nrows);
320 715 : std::vector<double> VT(ncols*ncols);
321 715 : std::vector<int> iwork(8*nsv);
322 :
323 : // This optimizes the size of the work array used in lapack singular value decomposition
324 715 : int lwork=-1;
325 715 : std::vector<double> work(1);
326 715 : plumed_lapack_dgesdd( "A", &nrows, &ncols, da.data(), &nrows, S.data(), U.data(), &nrows, VT.data(), &ncols, work.data(), &lwork, iwork.data(), &info );
327 : //if(info!=0) return info;
328 715 : if(info!=0) {
329 0 : log.printf("info:", info);
330 : }
331 :
332 : // Retrieve correct sizes for work and rellocate
333 715 : lwork=(int) work[0];
334 715 : work.resize(lwork);
335 :
336 : // This does the singular value decomposition
337 715 : plumed_lapack_dgesdd( "A", &nrows, &ncols, da.data(), &nrows, S.data(), U.data(), &nrows, VT.data(), &ncols, work.data(), &lwork, iwork.data(), &info );
338 : //if(info!=0) return info;
339 715 : if(info!=0) {
340 0 : log.printf("info:", info);
341 : }
342 :
343 :
344 : // get U and VT in form of 2D vector (U_, VT_)
345 715 : std::vector<std::vector<double>> U_(nrows, std::vector<double>(nrows,0));
346 715 : std::vector<std::vector<double>> VT_(ncols, std::vector<double>(ncols,0));
347 :
348 : int c=0;
349 :
350 2860 : for(int i=0; i<nrows; ++i) {
351 8580 : for(int j=0; j<nrows; ++j) {
352 6435 : U_[j][i] = U[c];
353 6435 : c += 1;
354 : }
355 : }
356 : c = 0; // note! its [j][i] not [i][j]
357 2860 : for(int i=0; i<ncols; ++i) {
358 8580 : for(int j=0; j<ncols; ++j) {
359 6435 : VT_[j][i] = VT[c];
360 6435 : c += 1;
361 : }
362 : }
363 : c=0; // note! its [j][i] not [i][j]
364 :
365 :
366 : // calculate determinants
367 715 : double det_u = determinant(nrows, &U_);
368 715 : double det_vt = determinant(ncols, &VT_);
369 :
370 : // check!
371 715 : if (det_u * det_vt < 0.0) {
372 1144 : for(int i=0; i<nrows; ++i) {
373 858 : U_[i][nrows-1] *= -1;
374 : }
375 : }
376 :
377 :
378 : //Matrix<double> rotation(3,3);
379 715 : rotation.resize(3,3);
380 : Matrix<double> u(3,3), vt(3,3);
381 2860 : for(int i=0; i<3; ++i) {
382 8580 : for(int j=0; j<3; ++j) {
383 6435 : u(i,j)=U_[i][j];
384 6435 : vt(i,j)=VT_[i][j];
385 : }
386 : }
387 :
388 : // get rotation matrix
389 715 : mult(u, vt, rotation);
390 :
391 1430 : }
392 :
393 :
394 : // calculates linear projection
395 715 : double position_linear_proj::cal_position_linear_proj() {
396 :
397 : unsigned N=getNumberOfAtoms();
398 :
399 : Matrix<double> rotated_obj(N,3);
400 : // rotate the object
401 715 : mult(mobile_str, rotation, rotated_obj);
402 :
403 : // compute the displacement
404 715 : std::vector<double> disp(N*3);
405 : unsigned c=0;
406 17160 : for(unsigned int i=0; i<N; ++i) {
407 65780 : for(int j=0; j<3; ++j) {
408 49335 : disp[c] = (rotated_obj(i,j)-ref_str(i,j));
409 49335 : c+=1;
410 : }
411 : }
412 :
413 : //double proj_val = dotProduct(disp, linear_coeffs);
414 : double proj_val = 0.0;
415 50050 : for(unsigned int i=0; i<N*3; ++i) {
416 49335 : proj_val += (linear_coeffs[i]*disp[i]);
417 : }
418 :
419 715 : return proj_val;
420 : }
421 :
422 : // numeric gradient
423 25 : void position_linear_proj::numeric_grad() {
424 : // This function performs numerical derivative.
425 : unsigned N=getNumberOfAtoms();
426 :
427 : unsigned stride;
428 : unsigned rank;
429 25 : if(serial) {
430 : // when using components the parallelisation do not work
431 : stride=1;
432 : rank=0;
433 : } else {
434 25 : stride=comm.Get_size();
435 25 : rank=comm.Get_rank();
436 : }
437 :
438 255 : for(unsigned i=rank; i<N; i+=stride) {
439 920 : for (unsigned j=0; j<3; ++j) {
440 :
441 690 : mobile_str(i,j) += delta;
442 690 : kabsch_rot_mat();
443 690 : derv_numeric(i,j) = ((cal_position_linear_proj() - proj)/delta);
444 :
445 690 : mobile_str(i,j) -= delta;
446 : }
447 :
448 : }
449 :
450 25 : if(!serial) {
451 25 : if(!derv_numeric.getVector().empty()) {
452 25 : comm.Sum(&derv_numeric(0,0), derv_numeric.getVector().size());
453 : }
454 : }
455 :
456 :
457 600 : for(unsigned i=0; i<N; ++i) {
458 575 : Vector vi(derv_numeric(i,0), derv_numeric(i,1), derv_numeric(i,2) );
459 575 : setAtomsDerivatives(i, vi);
460 : }
461 :
462 : // clear the matrix (very important step!!)
463 25 : derv_numeric *= 0;
464 25 : }
465 :
466 :
467 : // calculator
468 25 : void position_linear_proj::calculate() {
469 :
470 25 : if(pbc) {
471 25 : makeWhole();
472 : }
473 : unsigned N=getNumberOfAtoms();
474 :
475 : mobile_str.resize(N,3);
476 :
477 : // load the mobile str
478 600 : for(unsigned int i=0; i<N; ++i) {
479 575 : Vector pos=getPosition(i); // const PLMD::Vector
480 2300 : for(unsigned j=0; j<3; ++j) {
481 1725 : mobile_str(i,j) = pos[j];
482 : }
483 : }
484 :
485 : // translating the structure to center of geometry
486 25 : double center_of_geometry[3]= {0.0, 0.0, 0.0};
487 :
488 600 : for(unsigned int i=0; i<N; ++i) {
489 575 : center_of_geometry[0] += mobile_str(i,0);
490 575 : center_of_geometry[1] += mobile_str(i,1);
491 575 : center_of_geometry[2] += mobile_str(i,2);
492 : }
493 :
494 600 : for(unsigned int i=0; i<N; ++i) {
495 2300 : for(int j=0; j<3; ++j) {
496 1725 : mobile_str(i,j) -= (center_of_geometry[j]/N);
497 : }
498 : }
499 :
500 25 : kabsch_rot_mat();
501 25 : proj = cal_position_linear_proj();
502 :
503 25 : numeric_grad();
504 25 : setBoxDerivativesNoPbc();
505 25 : setValue(proj);
506 :
507 :
508 25 : }
509 :
510 : }
511 : }
512 :
513 :
514 :
|
