All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Modules Pages
Kearsley.cpp
Go to the documentation of this file.
1 /* +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
2  Copyright (c) 2013 The plumed team
3  (see the PEOPLE file at the root of the distribution for a list of names)
4 
5  See http://www.plumed-code.org for more information.
6 
7  This file is part of plumed, version 2.0.
8 
9  plumed is free software: you can redistribute it and/or modify
10  it under the terms of the GNU Lesser General Public License as published by
11  the Free Software Foundation, either version 3 of the License, or
12  (at your option) any later version.
13 
14  plumed is distributed in the hope that it will be useful,
15  but WITHOUT ANY WARRANTY; without even the implied warranty of
16  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
17  GNU Lesser General Public License for more details.
18 
19  You should have received a copy of the GNU Lesser General Public License
20  along with plumed. If not, see <http://www.gnu.org/licenses/>.
21 +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ */
22 #include "Kearsley.h"
23 #include <cmath>
24 #include <iostream>
25 #include <cstdlib>
26 #include "Matrix.h"
27 #include "Tensor.h"
28 #include "Log.h"
29 #include "Matrix.h"
30 #include "Random.h"
31 
32 using namespace std;
33 namespace PLMD{
34 
35 // put some notes
36 
37 Kearsley::Kearsley(const vector<Vector> &p0, const vector<Vector> &p1, const vector<double> &align, Log* &log):
38  log(log),
39  p0(p0),
40  p1(p1),
41  align(align),
42  com0_is_removed(false),
43  com1_is_removed(false),
44  err(0.0)
45 {
46  // now make an initial allocation
47 // int n=p0.size();
48  // eventually here one should make a "hard" resize of all the structures
49 // log->printf("Reallocating a size of %d atoms to kearsley structure\n",n);
50 
51 // try{
52 // diff.resize(n);
53 // }catch(bad_alloc&) {
54 // cerr<<"Cannot allocate the vector in Kearsley"<<endl;
55 // }
56 // try{
57 // diff.resize(n);
58 // }catch(bad_alloc&) {
59 // cerr<<"Cannot allocate the vector in Kearsley"<<endl;
60 // }
61 
62 }
63 // do the alignment
64 
65 double Kearsley::calculate(bool rmsd) {
66  // just an ad-hoc scaling factor
67  double myscale=1.;
68  // basic sanity check
69  if(p0.size()!=p1.size() || p1.size()!=align.size()){
70  cerr<<"Kearsley: looks like you have not properly allocated the vectors: the two frames have different size"<<endl;
71  cerr<<"size of p0 is :"<<p0.size()<<endl;
72  cerr<<"size of p1 is :"<<p1.size()<<endl;
73  cerr<<"size of align is :"<<align.size()<<endl;
74  exit(0);
75  }
76  if(p0.empty() || p1.empty() ){
77  cerr<<"Kearsley: looks like you have not properly allocated the vectors: they do not contain anything"<<endl;
78  exit(0);
79  }
80  Vector rr1,rr0;
81  Vector4d q;
82  double dddq[4][4][4],gamma[3][3][3],rrsq;
83  Matrix<double> m=Matrix<double>(4,4);
84 // double dm_r1[4][4][3],dm_r0[4][4][3];
85  Vector dm_r1[4][4];
86  Vector dm_r0[4][4];
87  Tensor pi1,pi0;
88  Tensor d;
89  Tensor dinv;
90  Vector xx;
91  double totalign=0., s, tmp1, err;
92  unsigned i,j,k,l,ii,ll,jj,mm,n,nn,iii;
93  bool verbose=false;
94 
95  vector<int> alignmap; // on the fly map done for optimization
96 
97  unsigned natoms=p0.size();
98 
99 
100  // calculate coms
101 
102  for(i=0;i<natoms;i++){
103  if (align[i]>0.){
104  alignmap.push_back(i);
105  totalign+=align[i];
106  }
107  if (align[i]<0.){cerr<<"FOUND ALIGNMENT WEIGHT NEGATIVE!"<<endl;exit(0);};
108 
109  }
110 
111 
112  // later will be implemented something for optimizing this piece of crap
113 
114  com0_is_removed=false;
115  com1_is_removed=false;
116 
117  bool do_center=true; // keep it for legacy code compatibility
118 
119  if(com0_is_removed==false){
120 
121  xx.zero();
122 
123  if(do_center) {// if you dont need to center no prob...
124  for(i=0;i<alignmap.size();i++){
125  xx+=p0[alignmap[i]]*align[alignmap[i]];
126  }
127  xx/=(totalign);
128  }
129 
130  com0=xx;
131 
132  if (p0reset.empty()){p0reset.resize(natoms);}
133  for(i=0;i<natoms;i++){
134  p0reset[i]=p0[i]-xx;
135  }
136  com0_is_removed=true;
137 
138  if (verbose){
139  log->printf("P0 RESET\n");
140  for(i=0;i<natoms;i++){
141  log->printf("ATOM %6d C ALA 1 %8.3f%8.3f%8.3f 1.00 1.00\n",i+1,p0reset[i][0]/myscale,p0reset[i][1]/myscale,p0reset[i][2]/myscale);
142  }
143  log->printf("END\n");
144  }
145 
146  }
147 
148  if(com1_is_removed==false){
149 
150  xx.zero();
151 
152  if(do_center) {// if you dont need to center no prob...
153  for(i=0;i<alignmap.size();i++){
154  xx+=p1[alignmap[i]]*align[alignmap[i]];
155  }
156  xx/=(totalign);
157  }
158 
159  com1=xx;
160 
161  if (p1reset.empty()){p1reset.resize(natoms);}
162 
163  for(i=0;i<natoms;i++){
164  p1reset[i]=p1[i]-xx;
165  }
166  com1_is_removed=true;
167 
168  if(verbose){
169  log->printf("P1 RESET\n");
170  for(i=0;i<natoms;i++){
171  log->printf("ATOM %6d C ALA 1 %8.3f%8.3f%8.3f 1.00 1.00\n",i+1,p1reset[i][0]/myscale,p1reset[i][1]/myscale,p1reset[i][2]/myscale);
172  }
173  log->printf("END\n");
174  }
175 
176  }
177 
178  bool fake=false;
179  if(fake){
180  // case of trivial alignment
181  // set rotmat
182  rotmat0on1=Tensor::identity();
183  if (p0reset.size()==0){p0reset.resize(natoms);}
184  if (p1reset.size()==0){p1reset.resize(natoms);}
185 
186  derrdp0.resize(natoms);
187  derrdp1.resize(natoms);
188  dmatdp0.resize(3*3*3*natoms);for(i=0;i<dmatdp0.size();i++)dmatdp0[i]=0.;
189  dmatdp1.resize(3*3*3*natoms);for(i=0;i<dmatdp1.size();i++)dmatdp1[i]=0.;
190 
191  err=0.;
192  for(i=0;i<natoms;i++){
193  if(align[i]>0.)err+=align[i]*modulo2(p0reset[i]-p1reset[i]);
194  }
195 
196  return 0.;
197  }
198  //
199  // CLEAN M MATRIX
200 
201  m=0.;
202 
203  // ASSIGN MATRIX ELEMENTS USING ONLY THE ATOMS INVOLVED IN ALIGNMENT
204 
205  for(i=0;i<alignmap.size();i++){
206 
207  k=alignmap[i];
208  tmp1=sqrt(align[k]/totalign);
209 
210  // adopt scaled coordinates
211 
212  rr1=p1reset[k]*tmp1;
213  rr0=p0reset[k]*tmp1;
214 
215  rrsq=modulo2(rr0)+modulo2(rr1);
216 
217  m[0][0] += rrsq+2.*(-rr0[0]*rr1[0]-rr0[1]*rr1[1]-rr0[2]*rr1[2]);
218  m[1][1] += rrsq+2.*(-rr0[0]*rr1[0]+rr0[1]*rr1[1]+rr0[2]*rr1[2]);
219  m[2][2] += rrsq+2.*(+rr0[0]*rr1[0]-rr0[1]*rr1[1]+rr0[2]*rr1[2]);
220  m[3][3] += rrsq+2.*(+rr0[0]*rr1[0]+rr0[1]*rr1[1]-rr0[2]*rr1[2]);
221  m[0][1] += 2.*(-rr0[1]*rr1[2]+rr0[2]*rr1[1]);
222  m[0][2] += 2.*( rr0[0]*rr1[2]-rr0[2]*rr1[0]);
223  m[0][3] += 2.*(-rr0[0]*rr1[1]+rr0[1]*rr1[0]);
224  m[1][2] -= 2.*( rr0[0]*rr1[1]+rr0[1]*rr1[0]);
225  m[1][3] -= 2.*( rr0[0]*rr1[2]+rr0[2]*rr1[0]);
226  m[2][3] -= 2.*( rr0[1]*rr1[2]+rr0[2]*rr1[1]);
227 
228  };
229  m[1][0] = m[0][1];
230  m[2][0] = m[0][2];
231  m[2][1] = m[1][2];
232  m[3][0] = m[0][3];
233  m[3][1] = m[1][3];
234  m[3][2] = m[2][3];
235 
236  // diagonalize the 4x4 matrix
237 
238  vector<double> eigenvals;
239  Matrix<double> eigenvecs;
240 
241  int diagerror=diagMat(m, eigenvals, eigenvecs );
242 
243  if (diagerror!=0){cerr<<"DIAGONALIZATION FAILED WITH ERROR CODE "<<diagerror<<endl;exit(0);}
244 
245  s=1.0;
246  if(eigenvecs(0,0)<0.)s=-1.;//correct for negative values (?)
247  // eigenvecs are in rows!!
248 
249  q[0]=s*eigenvecs(0,0);
250  q[1]=s*eigenvecs(0,1);
251  q[2]=s*eigenvecs(0,2);
252  q[3]=s*eigenvecs(0,3);
253  err=eigenvals[0];
254 
255  //log->printf(" ERR: %20.10f \n",err);
256 
257  if(verbose){
258  log->printf(" ERR: %f \n",err);
259  for (i=0;i<4;i++){
260  log->printf(" EIGENVALS: %f \n",eigenvals[i]);
261  }
262  }
263 
264  if(abs(eigenvals[0]-eigenvals[1])<1.e-8){
265  cerr<<"DIAGONALIZATION: NON UNIQUE SOLUTION"<<endl;exit(0);
266  }
267 
268  /*
269  * the ROTATION matrix
270  */
271 
272  d[0][0]=q[0]*q[0]+q[1]*q[1]-q[2]*q[2]-q[3]*q[3] ;
273  d[1][0]=2.0*(q[1]*q[2]-q[0]*q[3]);
274  d[2][0]=2.0*(q[1]*q[3]+q[0]*q[2]);
275  d[0][1]=2.0*(q[1]*q[2]+q[0]*q[3]);
276  d[1][1]=q[0]*q[0]+q[2]*q[2]-q[1]*q[1]-q[3]*q[3];
277  d[2][1]=2.0*(q[2]*q[3]-q[0]*q[1]);
278  d[0][2]=2.0*(q[1]*q[3]-q[0]*q[2]);
279  d[1][2]=2.0*(q[2]*q[3]+q[0]*q[1]);
280  d[2][2]=q[0]*q[0]+q[3]*q[3]-q[1]*q[1]-q[2]*q[2];
281 
282  /*
283  * first derivative in perturbation theory : derivative of the rotation matrix respect to the
284  * quternion vectors
285  */
286 
287  dddq[0][0][0]= 2.0*q[0];
288  dddq[1][0][0]=-2.0*q[3];
289  dddq[2][0][0]= 2.0*q[2];
290  dddq[0][1][0]= 2.0*q[3];
291  dddq[1][1][0]= 2.0*q[0];
292  dddq[2][1][0]=-2.0*q[1];
293  dddq[0][2][0]=-2.0*q[2];
294  dddq[1][2][0]= 2.0*q[1];
295  dddq[2][2][0]= 2.0*q[0];
296 
297  dddq[0][0][1]= 2.0*q[1];
298  dddq[1][0][1]= 2.0*q[2];
299  dddq[2][0][1]= 2.0*q[3];
300  dddq[0][1][1]= 2.0*q[2];
301  dddq[1][1][1]=-2.0*q[1];
302  dddq[2][1][1]=-2.0*q[0];
303  dddq[0][2][1]= 2.0*q[3];
304  dddq[1][2][1]= 2.0*q[0];
305  dddq[2][2][1]=-2.0*q[1];
306 
307  dddq[0][0][2]=-2.0*q[2];
308  dddq[1][0][2]= 2.0*q[1];
309  dddq[2][0][2]= 2.0*q[0];
310  dddq[0][1][2]= 2.0*q[1];
311  dddq[1][1][2]= 2.0*q[2];
312  dddq[2][1][2]= 2.0*q[3];
313  dddq[0][2][2]=-2.0*q[0];
314  dddq[1][2][2]= 2.0*q[3];
315  dddq[2][2][2]=-2.0*q[2];
316 
317  dddq[0][0][3]=-2.0*q[3];
318  dddq[1][0][3]=-2.0*q[0];
319  dddq[2][0][3]= 2.0*q[1];
320  dddq[0][1][3]= 2.0*q[0];
321  dddq[1][1][3]=-2.0*q[3];
322  dddq[2][1][3]= 2.0*q[2];
323  dddq[0][2][3]= 2.0*q[1];
324  dddq[1][2][3]= 2.0*q[2];
325  dddq[2][2][3]= 2.0*q[3];
326 
327  /*
328  * Build gamma 3x3x3 matrix
329  */
330  for(i=0;i<3;i++){ //direction
331  for(j=0;j<3;j++){ //direction
332  for(k=0;k<3;k++){ //eigenvector number
333  gamma[i][j][k]=0.0;
334  for(l=0;l<4;l++){ //components of each eigenvector in pert. series
335  if(abs(eigenvals[0]-eigenvals[k+1])<1.e-8){
336  log->printf(" FOUND DEGENERACY IN RMSD_ESS ROUTINE \n");
337  log->printf(" I'm DYING....\n");
338  log->printf(" COPYING STACK HERE \n");
339  log->printf(" P0\n");
340  for(ll=0;ll<natoms;ll++)log->printf(" %f %f %f \n",p0reset[ll][0],p0reset[ll][1],p0reset[ll][2]);
341  log->printf(" P1\n");
342  for(ll=0;ll<natoms;ll++)log->printf(" %f %f %f \n",p1reset[ll][0],p1reset[ll][1],p1reset[ll][2]);
343  exit(0);
344  }
345  else{
346  gamma[i][j][k] += dddq[i][j][l]*eigenvecs(k+1,l)/(eigenvals[0]-eigenvals[k+1]);
347  }
348  }
349  //log->printf("GAMMA %2d %2d %2d V %12.6f\n",i,j,k,gamma[i][j][k]);
350  }
351  }
352  }
353 
354  // allocate various arrays
355 
356  dmatdp1.resize(3*3*3*natoms);
357  for(i=0;i<dmatdp1.size();i++)dmatdp1[i]=0.;
358  dmatdp0.resize(3*3*3*natoms);
359  for(i=0;i<dmatdp0.size();i++)dmatdp0[i]=0.;
360 
361  vector<double> dd_dr_temp;dd_dr_temp.resize(natoms);
362 
363  vector<Vector> derr_dr1;
364  derr_dr1.resize(natoms);
365  vector<Vector> derr_dr0;
366  derr_dr0.resize(natoms);
367  vector<Vector> array_3_n;
368  array_3_n.resize(natoms);
369 
370 
371  /*
372  * Table of Derivative of the quaternion matrix respect to atom position: needed only if simple
373  * alignment is required and no correction respect to the rotation matrix is wanted
374  */
375 
376  for(iii=0;iii<alignmap.size();iii++){
377 
378  i=alignmap[iii];
379  tmp1=sqrt(align[i]/totalign);
380 
381  // once again: derivative respect to scaled distance
382 
383  rr1=2.*p1reset[i]*tmp1;
384  rr0=2.*p0reset[i]*tmp1;
385 
386 
387  dm_r1 [0][0][0]=(rr1[0]-rr0[0]);
388  dm_r1 [0][0][1]=(rr1[1]-rr0[1]);
389  dm_r1 [0][0][2]=(rr1[2]-rr0[2]);
390 
391  dm_r1 [0][1][0]=0.;
392  dm_r1 [0][1][1]= rr0[2];
393  dm_r1 [0][1][2]=-rr0[1];
394 
395  dm_r1 [0][2][0]=-rr0[2];
396  dm_r1 [0][2][1]= 0.;
397  dm_r1 [0][2][2]= rr0[0];
398 
399  dm_r1 [0][3][0]= rr0[1];
400  dm_r1 [0][3][1]=-rr0[0];
401  dm_r1 [0][3][2]= 0.;
402 
403  dm_r1 [1][1][0]=(rr1[0]-rr0[0]);
404  dm_r1 [1][1][1]=(rr1[1]+rr0[1]);
405  dm_r1 [1][1][2]=(rr1[2]+rr0[2]);
406 
407  dm_r1 [1][2][0]=-rr0[1];
408  dm_r1 [1][2][1]=-rr0[0];
409  dm_r1 [1][2][2]= 0.;
410 
411  dm_r1 [1][3][0]=-rr0[2];
412  dm_r1 [1][3][1]= 0.;
413  dm_r1 [1][3][2]=-rr0[0];
414 
415  dm_r1 [2][2][0]=(rr1[0]+rr0[0]);
416  dm_r1 [2][2][1]=(rr1[1]-rr0[1]);
417  dm_r1 [2][2][2]=(rr1[2]+rr0[2]);
418 
419  dm_r1 [2][3][0]=0.;
420  dm_r1 [2][3][1]=-rr0[2];
421  dm_r1 [2][3][2]=-rr0[1];
422 
423  dm_r1 [3][3][0]=(rr1[0]+rr0[0]);
424  dm_r1 [3][3][1]=(rr1[1]+rr0[1]);
425  dm_r1 [3][3][2]=(rr1[2]-rr0[2]);
426  /*
427  derivative respec to to the other vector
428  */
429  dm_r0 [0][0][0]=-(rr1[0]-rr0[0]);
430  dm_r0 [0][0][1]=-(rr1[1]-rr0[1]);
431  dm_r0 [0][0][2]=-(rr1[2]-rr0[2]);
432 
433  dm_r0 [0][1][0]=0. ;
434  dm_r0 [0][1][1]=-rr1[2];
435  dm_r0 [0][1][2]=rr1[1];
436 
437  dm_r0 [0][2][0]= rr1[2];
438  dm_r0 [0][2][1]= 0.;
439  dm_r0 [0][2][2]=-rr1[0];
440 
441  dm_r0 [0][3][0]=-rr1[1] ;
442  dm_r0 [0][3][1]= rr1[0];
443  dm_r0 [0][3][2]= 0.;
444 
445  dm_r0 [1][1][0]=-(rr1[0]-rr0[0]);
446  dm_r0 [1][1][1]=(rr1[1]+rr0[1]);
447  dm_r0 [1][1][2]=(rr1[2]+rr0[2]);
448 
449  dm_r0 [1][2][0]=-rr1[1];
450  dm_r0 [1][2][1]=-rr1[0];
451  dm_r0 [1][2][2]= 0.;
452 
453  dm_r0 [1][3][0]=-rr1[2];
454  dm_r0 [1][3][1]= 0.;
455  dm_r0 [1][3][2]=-rr1[0];
456 
457  dm_r0 [2][2][0]=(rr1[0]+rr0[0]);
458  dm_r0 [2][2][1]=-(rr1[1]-rr0[1]);
459  dm_r0 [2][2][2]=(rr1[2]+rr0[2]);
460 
461  dm_r0 [2][3][0]=0.;
462  dm_r0 [2][3][1]=-rr1[2];
463  dm_r0 [2][3][2]=-rr1[1];
464 
465  dm_r0 [3][3][0]=(rr1[0]+rr0[0]);
466  dm_r0 [3][3][1]=(rr1[1]+rr0[1]);
467  dm_r0 [3][3][2]=-(rr1[2]-rr0[2]);
468  /*
469  * write the diagonal
470  */
471 
472  dm_r1[1][0]=dm_r1[0][1];
473  dm_r1[2][0]=dm_r1[0][2];
474  dm_r1[3][0]=dm_r1[0][3];
475  dm_r1[2][1]=dm_r1[1][2];
476  dm_r1[3][1]=dm_r1[1][3];
477  dm_r1[3][2]=dm_r1[2][3];
478 
479  dm_r0[1][0]=dm_r0[0][1];
480  dm_r0[2][0]=dm_r0[0][2];
481  dm_r0[3][0]=dm_r0[0][3];
482  dm_r0[2][1]=dm_r0[1][2];
483  dm_r0[3][1]=dm_r0[1][3];
484  dm_r0[3][2]=dm_r0[2][3];
485 
486 
487  //log->printf("DMDR0 ALIGN %f AT %d VAL %f\n",align[i],i,dm_r0[0][0][j]);
488 
489 
490  /*
491  * pi matrix : coefficents in per theory
492  */
493 
494  pi0.zero();
495  pi1.zero();
496  derr_dr1[i].zero();
497  derr_dr0[i].zero();
498 
499  for(k=0;k<4;k++){
500  for(l=0;l<4;l++){
501  derr_dr1[i]+=(q[k]*q[l])*dm_r1[l][k];
502  derr_dr0[i]+=(q[k]*q[l])*dm_r0[l][k];
503  for(mm=0;mm<3;mm++)for(j=0;j<3;j++){
504  pi0[mm][j]+=eigenvecs(mm+1,k)*dm_r0[l][k][j]*q[l];
505  pi1[mm][j]+=eigenvecs(mm+1,k)*dm_r1[l][k][j]*q[l];
506  };
507  };
508  };
509 /*
510  derr_dr1[i]/=totalign;
511  derr_dr0[i]/=totalign;
512 
513 
514 */
515 
516  for(j=0;j<3;j++){
517  for (k=0;k<3;k++){
518  for(l=0;l<3;l++){
519  int ind=j*3*3*natoms+k*3*natoms+l*natoms+i;
520  dmatdp0[ind]=0.;
521  dmatdp1[ind]=0.;
522  for(ii=0;ii<3;ii++){
523  dmatdp1[ind]+=gamma[j][k][ii]*pi1[ii][l];
524  dmatdp0[ind]+=gamma[j][k][ii]*pi0[ii][l];
525  }
526 
527  }
528 
529  }
530 
531  }
532  }
533 
534 
535 
536 
537  // end of the calculation of the derivative of the rotation matrix
538 
539  /*
540  * Now correct for center of mass: only if needed
541  *
542  */
543 
544  bool comcorr_r1=true;
545 
546  if(comcorr_r1){
547  for(k=0;k<alignmap.size();k++){
548  i=alignmap[k];
549  tmp1=sqrt(align[i]/totalign);
550  array_3_n[i]=tmp1*derr_dr1[i];
551  if(do_center){
552  for(jj=0;jj<alignmap.size();jj++){
553  j=alignmap[jj];
554  array_3_n[i]-=tmp1*(align[j]/totalign)*derr_dr1[j];
555  }
556  }
557  }
558  for(k=0;k<alignmap.size();k++){
559  i=alignmap[k];
560  derr_dr1[i]=array_3_n[i];
561  }
562  }
563 
564 
565  bool do_comcorr_r0=true;
566  //
567  // correction for r0 frame
568  //
569  if(do_comcorr_r0){
570  for(k=0;k<alignmap.size();k++){
571  i=alignmap[k];
572  tmp1=sqrt(align[i]/totalign);
573  array_3_n[i]=tmp1*derr_dr0[i];
574  if(do_center){
575  for(jj=0;jj<alignmap.size();jj++){
576  j=alignmap[jj];
577  array_3_n[i]-=tmp1*(align[j]/totalign)*derr_dr0[j];
578  }
579  }
580  }
581  for(k=0;k<alignmap.size();k++){
582  i=alignmap[k];
583  derr_dr0[i]=array_3_n[i];
584  }
585  }
586 
587 
588  bool do_der_r1=true;
589  bool do_der_r0=true;
590  bool do_der_rotmat=true;
591 
592  if(do_der_r1 && do_der_rotmat){
593  for(i=0;i<3;i++){
594  for(j=0;j<3;j++){
595  for(k=0;k<3;k++){
596  for(ll=0;ll<alignmap.size();ll++){
597  l=alignmap[ll];
598  int ind=i*3*3*natoms+j*3*natoms+k*natoms+l;
599  tmp1=sqrt(align[l]/totalign);
600  dd_dr_temp[l]=tmp1*dmatdp1[ind];
601  if(do_center){
602  for(nn=0;nn<alignmap.size();nn++){
603  n=alignmap[nn];
604  dd_dr_temp[l]-=dmatdp1[ind-l+n]*tmp1*align[n]/totalign;
605  }
606  }
607 
608  }
609  for(ll=0;ll<alignmap.size();ll++){
610  l=alignmap[ll];
611  int ind=i*3*3*natoms+j*3*natoms+k*natoms+l;
612  dmatdp1[ind]=dd_dr_temp[l];
613  }
614 
615  }
616  }
617  }
618 
619  }
620 
621  if(do_der_r0 && do_der_rotmat){
622  for(i=0;i<3;i++){
623  for(j=0;j<3;j++){
624  for(k=0;k<3;k++){
625  for(ll=0;ll<alignmap.size();ll++){
626  l=alignmap[ll];
627  int ind=i*3*3*natoms+j*3*natoms+k*natoms+l;
628  tmp1=sqrt(align[l]/totalign);
629  dd_dr_temp[l]=tmp1*dmatdp0[ind];
630  if(do_center){
631  for(nn=0;nn<alignmap.size();nn++){
632  n=alignmap[nn];
633  dd_dr_temp[l]-=dmatdp0[ind-l+n]*tmp1*align[n]/totalign;
634  }
635  }
636  }
637  for(ll=0;ll<alignmap.size();ll++){
638  l=alignmap[ll];
639  int ind=i*3*3*natoms+j*3*natoms+k*natoms+l;
640  dmatdp0[ind]=dd_dr_temp[l];
641  }
642  }
643  }
644  }
645  }
646 
647 
648  bool do_p1rotated=true;
649  if (do_p1rotated){
650  // resize if not allocated
651 
652  if(p1.size()!=p1rotated.size())p1rotated.resize(p1.size());
653 
654 // exit(0);
655 
656  for(i=0;i<natoms;i++) p1rotated[i]=matmul(d,p1reset[i]);
657 
658  // reallocate difference vectors
659  if(p1.size()!=diff1on0.size())diff1on0.resize(p1.size());
660  for(i=0;i<natoms;i++){
661  diff1on0[i]=p1rotated[i]-p0reset[i];
662  }
663 
664  if(verbose){
665  log->printf("P1-RESET-AND-ROTATED\n");
666  for(i=0;i<natoms;i++){
667  log->printf("ATOM %6d C ALA 2 %8.3f%8.3f%8.3f 1.00 1.00\n",i+1,p1rotated[i][0]/myscale,p1rotated[i][1]/myscale,p1rotated[i][2]/myscale);
668  }
669  log->printf("END\n");
670  log->printf("P0-RESET\n");
671  for(i=0;i<natoms;i++){
672  log->printf("ATOM %6d C ALA 2 %8.3f%8.3f%8.3f 1.00 1.00\n",i+1,p0reset[i][0]/myscale,p0reset[i][1]/myscale,p0reset[i][2]/myscale);
673  }
674  log->printf("END\n");
675  }
676 
677  }
678 
679 
680 
681  dinv=inverse(d);
682 
683  bool do_p0rotated=true;
684  if (do_p0rotated){
685  if(p0.size()!=p0rotated.size())p0rotated.resize(p0.size());
686  for(i=0;i<natoms;i++) p0rotated[i]=matmul(dinv,p0reset[i]);
687  if(p1.size()!=diff0on1.size())diff0on1.resize(p1.size());
688  for(i=0;i<natoms;i++) diff0on1[i]=p0rotated[i]-p1reset[i];
689  if(verbose){
690  log->printf("P0-RESET AND INVERSE ROTATED\n");
691  for(i=0;i<natoms;i++){
692  log->printf("ATOM %6d C ALA 1 %8.3f%8.3f%8.3f 1.00 1.00\n",i+1,p0rotated[i][0]/myscale,p0rotated[i][1]/myscale,p0rotated[i][2]/myscale);
693  }
694  log->printf("END\n");
695  log->printf("P1-RESET\n");
696  for(i=0;i<natoms;i++){
697  log->printf("ATOM %6d C ALA 2 %8.3f%8.3f%8.3f 1.00 1.00\n",i+1,p1reset[i][0]/myscale,p1reset[i][1]/myscale,p1reset[i][2]/myscale);
698  }
699  log->printf("END\n");
700  }
701  }
702  // copy on the official vectors:
703  rotmat0on1=d;
704  rotmat1on0=dinv;
705  derrdp0.resize(natoms);
706  derrdp0=derr_dr0;
707  derrdp1.resize(natoms);
708  derrdp1=derr_dr1;
709 
710  // now rescale accordingly for rmsd instead of msd
711  if(rmsd){
712  err=sqrt(err);
713  double tmp=0.5/err;
714  for(ii=0;ii<alignmap.size();ii++){
715  i=alignmap[ii];
716  derrdp0[i]*=tmp;
717  derrdp1[i]*=tmp;
718  }
719 
720  }
721 
722  return err;
723 
724 }
725 void Kearsley::assignP1(const std::vector<Vector> & p1) {
726  this->p1=p1;
727  com1_is_removed=false;
728 }
729 void Kearsley::assignP0(const std::vector<Vector> & p0) {
730  this->p0=p0;
731  com0_is_removed=false;
732 }
733 
734 void Kearsley::assignAlign(const std::vector<double> & align) {
735  this->align=align;
736 }
737 
738 void Kearsley::finiteDifferenceInterface(bool rmsd){
739 log->printf("Entering rmsd finite difference test system for kearsley\n");
740 log->printf("-------------------------------------------\n");
741 log->printf("TEST1: derivative of the value (derr_dr0/derr_dr1)\n");
742 //// test 1
743 unsigned i,j,l,m;
744 double step=1.e-6,olderr,delta;
745 // messing up a bit with align weights
746 double delta1;
747 vector<double> align1;
748 align1.resize(p0.size());
749 Random rnd;
750 for (i=0;i<p0.size();i++){
751  // draw a random number
752  delta=rnd.RandU01();
753  delta1=rnd.RandU01();
754  if(delta>delta1){
755  //if(delta>0.3){
756  align1[i]=delta;
757  }else{align1[i]=0.;};
758  // log->printf("ALIGN %d IS %8.3f\n",i,align1[i]);
759 }
760 assignAlign(align1);
761 //// get initial value of the error and derivative of it
762 olderr=calculate(rmsd);
763 log->printf("INITIAL ERROR VALUE: %e\n",olderr);
764 // store the matrix
765 Tensor old_rotmat0on1=rotmat0on1;
766 
767 //// get initial value of the error and derivative of it
768 
769 log->printf("TESTING: derrdp1 \n");
770 for(unsigned j=0;j<3;j++){
771  for(unsigned i=0;i<derrdp1.size();i++){
772  // random displacement
773  delta=(rnd.RandU01()-0.5)*2*step;
774  p1[i][j]+=delta;
775  com1_is_removed=false; // this is required whenever the assignment is not done with the methods
776  com0_is_removed=false; // this is required whenever the assignment is not done with the methods
777  err=calculate(rmsd);
778  //log->printf("INITIAL ERROR VALUE: %e NEW ERROR %e DELTA %e ELEM %d %d \n",olderr,err,delta,i,j );
779  p1[i][j]-=delta;
780  switch(j){
781  case 0:
782  log->printf("TESTING: X %4u ANAL %18.9f NUMER %18.9f DELTA %18.9f ALIGN %6.2f\n",i,derrdp1[i][j],(err-olderr)/delta,derrdp1[i][j]-(err-olderr)/delta,align[i]);break;
783  case 1:
784  log->printf("TESTING: Y %4u ANAL %18.9f NUMER %18.9f DELTA %18.9f ALIGN %6.2f\n",i,derrdp1[i][j],(err-olderr)/delta,derrdp1[i][j]-(err-olderr)/delta,align[i]);break;
785  case 2:
786  log->printf("TESTING: Z %4u ANAL %18.9f NUMER %18.9f DELTA %18.9f ALIGN %6.2f\n",i,derrdp1[i][j],(err-olderr)/delta,derrdp1[i][j]-(err-olderr)/delta,align[i]);break;
787 
788  }
789  }
790 }
791 //exit(0);
792 log->printf("TESTING: derrdp0 \n");
793 for(unsigned j=0;j<3;j++){
794  for(unsigned i=0;i<derrdp0.size();i++){
795  // random displacement
796  delta=(rnd.RandU01()-0.5)*2*step;
797  p0[i][j]+=delta;
798  com0_is_removed=false; // this is required whenever the assignment is not done with the methods
799  com1_is_removed=false; // this is required whenever the assignment is not done with the methods
800 
801  err=calculate(rmsd);
802  p0[i][j]-=delta;
803  switch(j){
804  case 0:
805  log->printf("TESTING: X %4u ANAL %18.9f NUMER %18.9f DELTA %18.9f ALIGN %6.2f\n",i,derrdp0[i][j],(err-olderr)/delta,derrdp0[i][j]-(err-olderr)/delta,align[i]);break;
806  case 1:
807  log->printf("TESTING: Y %4u ANAL %18.9f NUMER %18.9f DELTA %18.9f ALIGN %6.2f\n",i,derrdp0[i][j],(err-olderr)/delta,derrdp0[i][j]-(err-olderr)/delta,align[i]);break;
808  case 2:
809  log->printf("TESTING: Z %4u ANAL %18.9f NUMER %18.9f DELTA %18.9f ALIGN %6.2f\n",i,derrdp0[i][j],(err-olderr)/delta,derrdp0[i][j]-(err-olderr)/delta,align[i]);break;
810  }
811  }
812 }
813 
814 log->printf("TESTING: dmatdp0 \n");
815 for(l=0;l<3;l++){
816  for(m=0;m<3;m++){
817  for(j=0;j<3;j++){
818  for(i=0;i<p0.size();i++){
819  // random displacement
820  delta=(rnd.RandU01()-0.5)*2*step;
821  p0[i][j]+=delta;
822  com0_is_removed=false;
823  com1_is_removed=false;
824  calculate(rmsd);
825  p0[i][j]-=delta;
826  int ind=l*3*3*p0.size()+m*3*p0.size()+j*p0.size()+i;
827  switch(j){
828  case 0:
829  log->printf("TESTING: DMATDP0 [ %u ][ %u ]: X %u ANAL %18.9f NUMER %18.9f DELTA %18.9f ALIGN %6.2f\n",l,m,i,dmatdp0[ind],(rotmat0on1[l][m]- old_rotmat0on1[l][m])/delta,dmatdp0[ind]-(rotmat0on1[l][m]- old_rotmat0on1[l][m])/delta,align[i]);break;
830 
831  case 1:
832  log->printf("TESTING: DMATDP0 [ %u ][ %u ]: Y %u ANAL %18.9f NUMER %18.9f DELTA %18.9f ALIGN %6.2f\n",l,m,i,dmatdp0[ind],(rotmat0on1[l][m]- old_rotmat0on1[l][m])/delta,dmatdp0[ind]-(rotmat0on1[l][m]- old_rotmat0on1[l][m])/delta,align[i]);break;
833 
834  case 2:
835  log->printf("TESTING: DMATDP0 [ %u ][ %u ]: Z %u ANAL %18.9f NUMER %18.9f DELTA %18.9f ALIGN %6.2f\n",l,m,i,dmatdp0[ind],(rotmat0on1[l][m]- old_rotmat0on1[l][m])/delta,dmatdp0[ind]-(rotmat0on1[l][m]- old_rotmat0on1[l][m])/delta,align[i]);break;
836 
837  }
838  }
839  }
840  }
841 }
842 log->printf("TESTING: dmatdp1 \n");
843 for(l=0;l<3;l++){
844  for(m=0;m<3;m++){
845  for(j=0;j<3;j++){
846  for(i=0;i<p1.size();i++){
847  // random displacement
848  delta=(rnd.RandU01()-0.5)*2*step;
849  p1[i][j]+=delta;
850  com0_is_removed=false;
851  com1_is_removed=false;
852  calculate(rmsd);
853  p1[i][j]-=delta;
854  int ind=l*3*3*p1.size()+m*3*p1.size()+j*p1.size()+i;
855  switch(j){
856 
857  case 0:
858  log->printf("TESTING: DMATDP1 [ %u ][ %u ]: X %u ANAL %18.9f NUMER %18.9f DELTA %18.9f ALIGN %6.2f\n",l,m,i,dmatdp1[ind],(rotmat0on1[l][m]- old_rotmat0on1[l][m])/delta,dmatdp1[ind]-(rotmat0on1[l][m]- old_rotmat0on1[l][m])/delta,align[i]);break;
859 
860  case 1:
861  log->printf("TESTING: DMATDP1 [ %u ][ %u ]: Y %u ANAL %18.9f NUMER %18.9f DELTA %18.9f ALIGN %6.2f\n",l,m,i,dmatdp1[ind],(rotmat0on1[l][m]- old_rotmat0on1[l][m])/delta,dmatdp1[ind]-(rotmat0on1[l][m]- old_rotmat0on1[l][m])/delta,align[i]);break;
862 
863  case 2:
864  log->printf("TESTING: DMATDP1 [ %u ][ %u ]: Z %u ANAL %18.9f NUMER %18.9f DELTA %18.9f ALIGN %6.2f\n",l,m,i,dmatdp1[ind],(rotmat0on1[l][m]- old_rotmat0on1[l][m])/delta,dmatdp1[ind]-(rotmat0on1[l][m]- old_rotmat0on1[l][m])/delta,align[i]);break;
865 
866  }
867  }
868  }
869  }
870 }
871 
872  exit(0);
873 }
874 }
PLMD::TensorGeneric::zero
void zero()
set it to zero
Definition: Tensor.h:198
PLMD::Kearsley::p1rotated
std::vector< Vector > p1rotated
position rotated: p1
Definition: Kearsley.h:210
PLMD::Kearsley::dmatdp1
std::vector< double > dmatdp1
Definition: Kearsley.h:220
PLMD::inverse
TensorGeneric< 3, 3 > inverse(const TensorGeneric< 3, 3 > &t)
Definition: Tensor.h:386
PLMD::Kearsley::com1_is_removed
bool com1_is_removed
Definition: Kearsley.h:189
PLMD::TensorGeneric
Class implementing fixed size matrices of doubles.
Definition: Tensor.h:70
PLMD::matmul
TensorGeneric< n, l > matmul(const TensorGeneric< n, m > &a, const TensorGeneric< m, l > &b)
Definition: Tensor.h:343
PLMD::VectorGeneric
Class implementing fixed size vectors of doubles.
Definition: Vector.h:74
PLMD::Kearsley::rotmat1on0
Tensor rotmat1on0
Definition: Kearsley.h:212
PLMD::Kearsley::err
double err
error: the distance between two frames (might be rmsd/msd. See below)
Definition: Kearsley.h:193
std
STL namespace.
PLMD::Log
Class containing the log stream.
Definition: Log.h:35
n
void const char const char int * n
Definition: Matrix.h:42
PLMD::Kearsley::p1reset
std::vector< Vector > p1reset
position resetted wrt coms p1
Definition: Kearsley.h:206
PLMD::VectorGeneric::zero
void zero()
set it to zero
Definition: Vector.h:173
PLMD::Kearsley::com0_is_removed
bool com0_is_removed
Definition: Kearsley.h:188
PLMD::Kearsley::derrdp1
std::vector< Vector > derrdp1
derivatives: derivative of the error respect p1
Definition: Kearsley.h:216
PLMD::Kearsley::p0rotated
std::vector< Vector > p0rotated
position rotated: p0
Definition: Kearsley.h:208
PLMD::TensorGeneric::identity
static TensorGeneric< n, n > identity()
return an identity tensor
Definition: Tensor.h:318
PLMD::Kearsley::log
Log * log
general log reference that needs to be initialized when constructed
Definition: Kearsley.h:180
PLMD::Kearsley::assignP0
void assignP0(const std::vector< Vector > &p0)
switch the assignment of the structure p0 (e.g. at each md step)
Definition: Kearsley.cpp:729
PLMD::Kearsley::calculate
double calculate(bool rmsd)
Definition: Kearsley.cpp:65
PLMD::Kearsley::diff1on0
std::vector< Vector > diff1on0
displacement: the vector that goes from the p1 onto p0 (via inverse rotation)
Definition: Kearsley.h:197
PLMD::Kearsley::p0
std::vector< Vector > p0
position of atoms (first frame. In md is the running frame)
Definition: Kearsley.h:182
PLMD::Kearsley::align
std::vector< double > align
alignment weight: the rmsd/msd that it provides is only based on this scalar
Definition: Kearsley.h:186
PLMD::OFile::printf
int printf(const char *fmt,...)
Formatted output with explicit format - a la printf.
Definition: OFile.cpp:82
PLMD::delta
VectorGeneric< n > delta(const VectorGeneric< n > &v1, const VectorGeneric< n > &v2)
Definition: Vector.h:262
PLMD::Kearsley::p0reset
std::vector< Vector > p0reset
position resetted wrt coms p0
Definition: Kearsley.h:204
PLMD::Kearsley::com0
Vector com0
center of mass of p0
Definition: Kearsley.h:200
PLMD::Kearsley::assignP1
void assignP1(const std::vector< Vector > &p1)
derivatives: derivative of the error respect p1
Definition: Kearsley.cpp:725
PLMD::diagMat
int diagMat(const Matrix< T > &A, std::vector< double > &eigenvals, Matrix< double > &eigenvecs)
Definition: Matrix.h:203
PLMD::modulo2
double modulo2(const VectorGeneric< n > &v)
Definition: Vector.h:330
m
void const char const char int double int double double int int double int * m
Definition: Matrix.h:42
PLMD::Random::RandU01
double RandU01()
Definition: Random.cpp:56
PLMD::Kearsley::diff0on1
std::vector< Vector > diff0on1
displacement: the vector that goes from the p0 onto p1
Definition: Kearsley.h:195
PLMD::Kearsley::com1
Vector com1
center of mass of p1
Definition: Kearsley.h:202
PLMD::Kearsley::assignAlign
void assignAlign(const std::vector< double > &align)
transfer the alignment vector
Definition: Kearsley.cpp:734
PLMD::Kearsley::dmatdp0
std::vector< double > dmatdp0
derivative of the rotation matrix note the dimension 3x3 x 3 x N
Definition: Kearsley.h:219
PLMD::Kearsley::rotmat0on1
Tensor rotmat0on1
rotation matrices p0 on p1 and reverse (p1 over p0)
Definition: Kearsley.h:212
PLMD::Kearsley::p1
std::vector< Vector > p1
position of atoms (second frame. In md is the reference frame)
Definition: Kearsley.h:184
PLMD::Kearsley::derrdp0
std::vector< Vector > derrdp0
derivatives: derivative of the error respect p0
Definition: Kearsley.h:214
PLMD::Kearsley::finiteDifferenceInterface
void finiteDifferenceInterface(bool rmsd)
finite differences of all the relevant quantities: takes a bool which decides if giving back rmsd or ...
Definition: Kearsley.cpp:738
Hosted by GitHub   GitHub Logo doxygen 1.8.8