Line data Source code
1 : /* +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
2 : Copyright (c) 2011-2018 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.org for more information.
6 :
7 : This file is part of plumed, version 2.
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 0 : 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 0 : 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 0 : }
63 : // do the alignment
64 :
65 0 : double Kearsley::calculate(bool rmsd) {
66 : // just an ad-hoc scaling factor
67 0 : double myscale=1.;
68 : // basic sanity check
69 0 : if(p0.size()!=p1.size() || p1.size()!=align.size()) {
70 0 : cerr<<"Kearsley: looks like you have not properly allocated the vectors: the two frames have different size"<<endl;
71 0 : cerr<<"size of p0 is :"<<p0.size()<<endl;
72 0 : cerr<<"size of p1 is :"<<p1.size()<<endl;
73 0 : cerr<<"size of align is :"<<align.size()<<endl;
74 0 : exit(0);
75 : }
76 0 : if(p0.empty() || p1.empty() ) {
77 0 : cerr<<"Kearsley: looks like you have not properly allocated the vectors: they do not contain anything"<<endl;
78 0 : exit(0);
79 : }
80 0 : Vector rr1,rr0;
81 0 : Vector4d q;
82 : double dddq[4][4][4],gamma[3][3][3],rrsq;
83 0 : Matrix<double> m=Matrix<double>(4,4);
84 : // double dm_r1[4][4][3],dm_r0[4][4][3];
85 0 : Vector dm_r1[4][4];
86 0 : Vector dm_r0[4][4];
87 0 : Tensor pi1,pi0;
88 0 : Tensor d;
89 0 : Tensor dinv;
90 0 : Vector xx;
91 0 : double totalign=0., s, tmp1, err;
92 : unsigned i,j,k,l,ii,ll,jj,mm,n,nn,iii;
93 0 : bool verbose=false;
94 :
95 0 : vector<int> alignmap; // on the fly map done for optimization
96 :
97 0 : unsigned natoms=p0.size();
98 :
99 :
100 : // calculate coms
101 :
102 0 : for(i=0; i<natoms; i++) {
103 0 : if (align[i]>0.) {
104 0 : alignmap.push_back(i);
105 0 : totalign+=align[i];
106 : }
107 0 : 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 0 : com0_is_removed=false;
115 0 : com1_is_removed=false;
116 :
117 0 : bool do_center=true; // keep it for legacy code compatibility
118 :
119 0 : if(com0_is_removed==false) {
120 :
121 0 : xx.zero();
122 :
123 0 : if(do_center) {// if you dont need to center no prob...
124 0 : for(i=0; i<alignmap.size(); i++) {
125 0 : xx+=p0[alignmap[i]]*align[alignmap[i]];
126 : }
127 0 : xx/=(totalign);
128 : }
129 :
130 0 : com0=xx;
131 :
132 0 : if (p0reset.empty()) {p0reset.resize(natoms);}
133 0 : for(i=0; i<natoms; i++) {
134 0 : p0reset[i]=p0[i]-xx;
135 : }
136 0 : com0_is_removed=true;
137 :
138 0 : if (verbose) {
139 0 : log->printf("P0 RESET\n");
140 0 : for(i=0; i<natoms; i++) {
141 0 : log->printf("ATOM %6u 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 0 : log->printf("END\n");
144 : }
145 :
146 : }
147 :
148 0 : if(com1_is_removed==false) {
149 :
150 0 : xx.zero();
151 :
152 0 : if(do_center) {// if you dont need to center no prob...
153 0 : for(i=0; i<alignmap.size(); i++) {
154 0 : xx+=p1[alignmap[i]]*align[alignmap[i]];
155 : }
156 0 : xx/=(totalign);
157 : }
158 :
159 0 : com1=xx;
160 :
161 0 : if (p1reset.empty()) {p1reset.resize(natoms);}
162 :
163 0 : for(i=0; i<natoms; i++) {
164 0 : p1reset[i]=p1[i]-xx;
165 : }
166 0 : com1_is_removed=true;
167 :
168 0 : if(verbose) {
169 0 : log->printf("P1 RESET\n");
170 0 : for(i=0; i<natoms; i++) {
171 0 : log->printf("ATOM %6u 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 0 : log->printf("END\n");
174 : }
175 :
176 : }
177 :
178 0 : bool fake=false;
179 0 : if(fake) {
180 : // case of trivial alignment
181 : // set rotmat
182 0 : rotmat0on1=Tensor::identity();
183 0 : if (p0reset.size()==0) {p0reset.resize(natoms);}
184 0 : if (p1reset.size()==0) {p1reset.resize(natoms);}
185 :
186 0 : derrdp0.resize(natoms);
187 0 : derrdp1.resize(natoms);
188 0 : dmatdp0.resize(3*3*3*natoms); for(i=0; i<dmatdp0.size(); i++)dmatdp0[i]=0.;
189 0 : dmatdp1.resize(3*3*3*natoms); for(i=0; i<dmatdp1.size(); i++)dmatdp1[i]=0.;
190 :
191 0 : err=0.;
192 0 : for(i=0; i<natoms; i++) {
193 0 : if(align[i]>0.)err+=align[i]*modulo2(p0reset[i]-p1reset[i]);
194 : }
195 :
196 0 : return 0.;
197 : }
198 : //
199 : // CLEAN M MATRIX
200 :
201 0 : m=0.;
202 :
203 : // ASSIGN MATRIX ELEMENTS USING ONLY THE ATOMS INVOLVED IN ALIGNMENT
204 :
205 0 : for(i=0; i<alignmap.size(); i++) {
206 :
207 0 : k=alignmap[i];
208 0 : tmp1=sqrt(align[k]/totalign);
209 :
210 : // adopt scaled coordinates
211 :
212 0 : rr1=p1reset[k]*tmp1;
213 0 : rr0=p0reset[k]*tmp1;
214 :
215 0 : rrsq=modulo2(rr0)+modulo2(rr1);
216 :
217 0 : m[0][0] += rrsq+2.*(-rr0[0]*rr1[0]-rr0[1]*rr1[1]-rr0[2]*rr1[2]);
218 0 : m[1][1] += rrsq+2.*(-rr0[0]*rr1[0]+rr0[1]*rr1[1]+rr0[2]*rr1[2]);
219 0 : m[2][2] += rrsq+2.*(+rr0[0]*rr1[0]-rr0[1]*rr1[1]+rr0[2]*rr1[2]);
220 0 : m[3][3] += rrsq+2.*(+rr0[0]*rr1[0]+rr0[1]*rr1[1]-rr0[2]*rr1[2]);
221 0 : m[0][1] += 2.*(-rr0[1]*rr1[2]+rr0[2]*rr1[1]);
222 0 : m[0][2] += 2.*( rr0[0]*rr1[2]-rr0[2]*rr1[0]);
223 0 : m[0][3] += 2.*(-rr0[0]*rr1[1]+rr0[1]*rr1[0]);
224 0 : m[1][2] -= 2.*( rr0[0]*rr1[1]+rr0[1]*rr1[0]);
225 0 : m[1][3] -= 2.*( rr0[0]*rr1[2]+rr0[2]*rr1[0]);
226 0 : m[2][3] -= 2.*( rr0[1]*rr1[2]+rr0[2]*rr1[1]);
227 :
228 : };
229 0 : m[1][0] = m[0][1];
230 0 : m[2][0] = m[0][2];
231 0 : m[2][1] = m[1][2];
232 0 : m[3][0] = m[0][3];
233 0 : m[3][1] = m[1][3];
234 0 : m[3][2] = m[2][3];
235 :
236 : // diagonalize the 4x4 matrix
237 :
238 0 : vector<double> eigenvals;
239 0 : Matrix<double> eigenvecs;
240 :
241 0 : int diagerror=diagMat(m, eigenvals, eigenvecs );
242 :
243 0 : if (diagerror!=0) {cerr<<"DIAGONALIZATION FAILED WITH ERROR CODE "<<diagerror<<endl; exit(0);}
244 :
245 0 : s=1.0;
246 0 : if(eigenvecs(0,0)<0.)s=-1.;//correct for negative values (?)
247 : // eigenvecs are in rows!!
248 :
249 0 : q[0]=s*eigenvecs(0,0);
250 0 : q[1]=s*eigenvecs(0,1);
251 0 : q[2]=s*eigenvecs(0,2);
252 0 : q[3]=s*eigenvecs(0,3);
253 0 : err=eigenvals[0];
254 :
255 : //log->printf(" ERR: %20.10f \n",err);
256 :
257 0 : if(verbose) {
258 0 : log->printf(" ERR: %f \n",err);
259 0 : for (i=0; i<4; i++) {
260 0 : log->printf(" EIGENVALS: %f \n",eigenvals[i]);
261 : }
262 : }
263 :
264 0 : if(abs(eigenvals[0]-eigenvals[1])<1.e-8) {
265 0 : cerr<<"DIAGONALIZATION: NON UNIQUE SOLUTION"<<endl; exit(0);
266 : }
267 :
268 : /*
269 : * the ROTATION matrix
270 : */
271 :
272 0 : d[0][0]=q[0]*q[0]+q[1]*q[1]-q[2]*q[2]-q[3]*q[3] ;
273 0 : d[1][0]=2.0*(q[1]*q[2]-q[0]*q[3]);
274 0 : d[2][0]=2.0*(q[1]*q[3]+q[0]*q[2]);
275 0 : d[0][1]=2.0*(q[1]*q[2]+q[0]*q[3]);
276 0 : d[1][1]=q[0]*q[0]+q[2]*q[2]-q[1]*q[1]-q[3]*q[3];
277 0 : d[2][1]=2.0*(q[2]*q[3]-q[0]*q[1]);
278 0 : d[0][2]=2.0*(q[1]*q[3]-q[0]*q[2]);
279 0 : d[1][2]=2.0*(q[2]*q[3]+q[0]*q[1]);
280 0 : 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 0 : dddq[0][0][0]= 2.0*q[0];
288 0 : dddq[1][0][0]=-2.0*q[3];
289 0 : dddq[2][0][0]= 2.0*q[2];
290 0 : dddq[0][1][0]= 2.0*q[3];
291 0 : dddq[1][1][0]= 2.0*q[0];
292 0 : dddq[2][1][0]=-2.0*q[1];
293 0 : dddq[0][2][0]=-2.0*q[2];
294 0 : dddq[1][2][0]= 2.0*q[1];
295 0 : dddq[2][2][0]= 2.0*q[0];
296 :
297 0 : dddq[0][0][1]= 2.0*q[1];
298 0 : dddq[1][0][1]= 2.0*q[2];
299 0 : dddq[2][0][1]= 2.0*q[3];
300 0 : dddq[0][1][1]= 2.0*q[2];
301 0 : dddq[1][1][1]=-2.0*q[1];
302 0 : dddq[2][1][1]=-2.0*q[0];
303 0 : dddq[0][2][1]= 2.0*q[3];
304 0 : dddq[1][2][1]= 2.0*q[0];
305 0 : dddq[2][2][1]=-2.0*q[1];
306 :
307 0 : dddq[0][0][2]=-2.0*q[2];
308 0 : dddq[1][0][2]= 2.0*q[1];
309 0 : dddq[2][0][2]= 2.0*q[0];
310 0 : dddq[0][1][2]= 2.0*q[1];
311 0 : dddq[1][1][2]= 2.0*q[2];
312 0 : dddq[2][1][2]= 2.0*q[3];
313 0 : dddq[0][2][2]=-2.0*q[0];
314 0 : dddq[1][2][2]= 2.0*q[3];
315 0 : dddq[2][2][2]=-2.0*q[2];
316 :
317 0 : dddq[0][0][3]=-2.0*q[3];
318 0 : dddq[1][0][3]=-2.0*q[0];
319 0 : dddq[2][0][3]= 2.0*q[1];
320 0 : dddq[0][1][3]= 2.0*q[0];
321 0 : dddq[1][1][3]=-2.0*q[3];
322 0 : dddq[2][1][3]= 2.0*q[2];
323 0 : dddq[0][2][3]= 2.0*q[1];
324 0 : dddq[1][2][3]= 2.0*q[2];
325 0 : dddq[2][2][3]= 2.0*q[3];
326 :
327 : /*
328 : * Build gamma 3x3x3 matrix
329 : */
330 0 : for(i=0; i<3; i++) { //direction
331 0 : for(j=0; j<3; j++) { //direction
332 0 : for(k=0; k<3; k++) { //eigenvector number
333 0 : gamma[i][j][k]=0.0;
334 0 : for(l=0; l<4; l++) { //components of each eigenvector in pert. series
335 0 : if(abs(eigenvals[0]-eigenvals[k+1])<1.e-8) {
336 0 : log->printf(" FOUND DEGENERACY IN RMSD_ESS ROUTINE \n");
337 0 : log->printf(" I'm DYING....\n");
338 0 : log->printf(" COPYING STACK HERE \n");
339 0 : log->printf(" P0\n");
340 0 : for(ll=0; ll<natoms; ll++)log->printf(" %f %f %f \n",p0reset[ll][0],p0reset[ll][1],p0reset[ll][2]);
341 0 : log->printf(" P1\n");
342 0 : for(ll=0; ll<natoms; ll++)log->printf(" %f %f %f \n",p1reset[ll][0],p1reset[ll][1],p1reset[ll][2]);
343 0 : exit(0);
344 : }
345 : else {
346 0 : 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 0 : dmatdp1.resize(3*3*3*natoms);
357 0 : for(i=0; i<dmatdp1.size(); i++)dmatdp1[i]=0.;
358 0 : dmatdp0.resize(3*3*3*natoms);
359 0 : for(i=0; i<dmatdp0.size(); i++)dmatdp0[i]=0.;
360 :
361 0 : vector<double> dd_dr_temp; dd_dr_temp.resize(natoms);
362 :
363 0 : vector<Vector> derr_dr1;
364 0 : derr_dr1.resize(natoms);
365 0 : vector<Vector> derr_dr0;
366 0 : derr_dr0.resize(natoms);
367 0 : vector<Vector> array_3_n;
368 0 : 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 0 : for(iii=0; iii<alignmap.size(); iii++) {
377 :
378 0 : i=alignmap[iii];
379 0 : tmp1=sqrt(align[i]/totalign);
380 :
381 : // once again: derivative respect to scaled distance
382 :
383 0 : rr1=2.*p1reset[i]*tmp1;
384 0 : rr0=2.*p0reset[i]*tmp1;
385 :
386 :
387 0 : dm_r1 [0][0][0]=(rr1[0]-rr0[0]);
388 0 : dm_r1 [0][0][1]=(rr1[1]-rr0[1]);
389 0 : dm_r1 [0][0][2]=(rr1[2]-rr0[2]);
390 :
391 0 : dm_r1 [0][1][0]=0.;
392 0 : dm_r1 [0][1][1]= rr0[2];
393 0 : dm_r1 [0][1][2]=-rr0[1];
394 :
395 0 : dm_r1 [0][2][0]=-rr0[2];
396 0 : dm_r1 [0][2][1]= 0.;
397 0 : dm_r1 [0][2][2]= rr0[0];
398 :
399 0 : dm_r1 [0][3][0]= rr0[1];
400 0 : dm_r1 [0][3][1]=-rr0[0];
401 0 : dm_r1 [0][3][2]= 0.;
402 :
403 0 : dm_r1 [1][1][0]=(rr1[0]-rr0[0]);
404 0 : dm_r1 [1][1][1]=(rr1[1]+rr0[1]);
405 0 : dm_r1 [1][1][2]=(rr1[2]+rr0[2]);
406 :
407 0 : dm_r1 [1][2][0]=-rr0[1];
408 0 : dm_r1 [1][2][1]=-rr0[0];
409 0 : dm_r1 [1][2][2]= 0.;
410 :
411 0 : dm_r1 [1][3][0]=-rr0[2];
412 0 : dm_r1 [1][3][1]= 0.;
413 0 : dm_r1 [1][3][2]=-rr0[0];
414 :
415 0 : dm_r1 [2][2][0]=(rr1[0]+rr0[0]);
416 0 : dm_r1 [2][2][1]=(rr1[1]-rr0[1]);
417 0 : dm_r1 [2][2][2]=(rr1[2]+rr0[2]);
418 :
419 0 : dm_r1 [2][3][0]=0.;
420 0 : dm_r1 [2][3][1]=-rr0[2];
421 0 : dm_r1 [2][3][2]=-rr0[1];
422 :
423 0 : dm_r1 [3][3][0]=(rr1[0]+rr0[0]);
424 0 : dm_r1 [3][3][1]=(rr1[1]+rr0[1]);
425 0 : dm_r1 [3][3][2]=(rr1[2]-rr0[2]);
426 : /*
427 : derivative respec to to the other vector
428 : */
429 0 : dm_r0 [0][0][0]=-(rr1[0]-rr0[0]);
430 0 : dm_r0 [0][0][1]=-(rr1[1]-rr0[1]);
431 0 : dm_r0 [0][0][2]=-(rr1[2]-rr0[2]);
432 :
433 0 : dm_r0 [0][1][0]=0. ;
434 0 : dm_r0 [0][1][1]=-rr1[2];
435 0 : dm_r0 [0][1][2]=rr1[1];
436 :
437 0 : dm_r0 [0][2][0]= rr1[2];
438 0 : dm_r0 [0][2][1]= 0.;
439 0 : dm_r0 [0][2][2]=-rr1[0];
440 :
441 0 : dm_r0 [0][3][0]=-rr1[1] ;
442 0 : dm_r0 [0][3][1]= rr1[0];
443 0 : dm_r0 [0][3][2]= 0.;
444 :
445 0 : dm_r0 [1][1][0]=-(rr1[0]-rr0[0]);
446 0 : dm_r0 [1][1][1]=(rr1[1]+rr0[1]);
447 0 : dm_r0 [1][1][2]=(rr1[2]+rr0[2]);
448 :
449 0 : dm_r0 [1][2][0]=-rr1[1];
450 0 : dm_r0 [1][2][1]=-rr1[0];
451 0 : dm_r0 [1][2][2]= 0.;
452 :
453 0 : dm_r0 [1][3][0]=-rr1[2];
454 0 : dm_r0 [1][3][1]= 0.;
455 0 : dm_r0 [1][3][2]=-rr1[0];
456 :
457 0 : dm_r0 [2][2][0]=(rr1[0]+rr0[0]);
458 0 : dm_r0 [2][2][1]=-(rr1[1]-rr0[1]);
459 0 : dm_r0 [2][2][2]=(rr1[2]+rr0[2]);
460 :
461 0 : dm_r0 [2][3][0]=0.;
462 0 : dm_r0 [2][3][1]=-rr1[2];
463 0 : dm_r0 [2][3][2]=-rr1[1];
464 :
465 0 : dm_r0 [3][3][0]=(rr1[0]+rr0[0]);
466 0 : dm_r0 [3][3][1]=(rr1[1]+rr0[1]);
467 0 : dm_r0 [3][3][2]=-(rr1[2]-rr0[2]);
468 : /*
469 : * write the diagonal
470 : */
471 :
472 0 : dm_r1[1][0]=dm_r1[0][1];
473 0 : dm_r1[2][0]=dm_r1[0][2];
474 0 : dm_r1[3][0]=dm_r1[0][3];
475 0 : dm_r1[2][1]=dm_r1[1][2];
476 0 : dm_r1[3][1]=dm_r1[1][3];
477 0 : dm_r1[3][2]=dm_r1[2][3];
478 :
479 0 : dm_r0[1][0]=dm_r0[0][1];
480 0 : dm_r0[2][0]=dm_r0[0][2];
481 0 : dm_r0[3][0]=dm_r0[0][3];
482 0 : dm_r0[2][1]=dm_r0[1][2];
483 0 : dm_r0[3][1]=dm_r0[1][3];
484 0 : 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 0 : pi0.zero();
495 0 : pi1.zero();
496 0 : derr_dr1[i].zero();
497 0 : derr_dr0[i].zero();
498 :
499 0 : for(k=0; k<4; k++) {
500 0 : for(l=0; l<4; l++) {
501 0 : derr_dr1[i]+=(q[k]*q[l])*dm_r1[l][k];
502 0 : derr_dr0[i]+=(q[k]*q[l])*dm_r0[l][k];
503 0 : for(mm=0; mm<3; mm++)for(j=0; j<3; j++) {
504 0 : pi0[mm][j]+=eigenvecs(mm+1,k)*dm_r0[l][k][j]*q[l];
505 0 : 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 0 : for(j=0; j<3; j++) {
517 0 : for (k=0; k<3; k++) {
518 0 : for(l=0; l<3; l++) {
519 0 : int ind=j*3*3*natoms+k*3*natoms+l*natoms+i;
520 0 : dmatdp0[ind]=0.;
521 0 : dmatdp1[ind]=0.;
522 0 : for(ii=0; ii<3; ii++) {
523 0 : dmatdp1[ind]+=gamma[j][k][ii]*pi1[ii][l];
524 0 : 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 0 : bool comcorr_r1=true;
545 :
546 0 : if(comcorr_r1) {
547 0 : for(k=0; k<alignmap.size(); k++) {
548 0 : i=alignmap[k];
549 0 : tmp1=sqrt(align[i]/totalign);
550 0 : array_3_n[i]=tmp1*derr_dr1[i];
551 0 : if(do_center) {
552 0 : for(jj=0; jj<alignmap.size(); jj++) {
553 0 : j=alignmap[jj];
554 0 : array_3_n[i]-=tmp1*(align[j]/totalign)*derr_dr1[j];
555 : }
556 : }
557 : }
558 0 : for(k=0; k<alignmap.size(); k++) {
559 0 : i=alignmap[k];
560 0 : derr_dr1[i]=array_3_n[i];
561 : }
562 : }
563 :
564 :
565 0 : bool do_comcorr_r0=true;
566 : //
567 : // correction for r0 frame
568 : //
569 0 : if(do_comcorr_r0) {
570 0 : for(k=0; k<alignmap.size(); k++) {
571 0 : i=alignmap[k];
572 0 : tmp1=sqrt(align[i]/totalign);
573 0 : array_3_n[i]=tmp1*derr_dr0[i];
574 0 : if(do_center) {
575 0 : for(jj=0; jj<alignmap.size(); jj++) {
576 0 : j=alignmap[jj];
577 0 : array_3_n[i]-=tmp1*(align[j]/totalign)*derr_dr0[j];
578 : }
579 : }
580 : }
581 0 : for(k=0; k<alignmap.size(); k++) {
582 0 : i=alignmap[k];
583 0 : derr_dr0[i]=array_3_n[i];
584 : }
585 : }
586 :
587 :
588 0 : bool do_der_r1=true;
589 0 : bool do_der_r0=true;
590 0 : bool do_der_rotmat=true;
591 :
592 0 : if(do_der_r1 && do_der_rotmat) {
593 0 : for(i=0; i<3; i++) {
594 0 : for(j=0; j<3; j++) {
595 0 : for(k=0; k<3; k++) {
596 0 : for(ll=0; ll<alignmap.size(); ll++) {
597 0 : l=alignmap[ll];
598 0 : int ind=i*3*3*natoms+j*3*natoms+k*natoms+l;
599 0 : tmp1=sqrt(align[l]/totalign);
600 0 : dd_dr_temp[l]=tmp1*dmatdp1[ind];
601 0 : if(do_center) {
602 0 : for(nn=0; nn<alignmap.size(); nn++) {
603 0 : n=alignmap[nn];
604 0 : dd_dr_temp[l]-=dmatdp1[ind-l+n]*tmp1*align[n]/totalign;
605 : }
606 : }
607 :
608 : }
609 0 : for(ll=0; ll<alignmap.size(); ll++) {
610 0 : l=alignmap[ll];
611 0 : int ind=i*3*3*natoms+j*3*natoms+k*natoms+l;
612 0 : dmatdp1[ind]=dd_dr_temp[l];
613 : }
614 :
615 : }
616 : }
617 : }
618 :
619 : }
620 :
621 0 : if(do_der_r0 && do_der_rotmat) {
622 0 : for(i=0; i<3; i++) {
623 0 : for(j=0; j<3; j++) {
624 0 : for(k=0; k<3; k++) {
625 0 : for(ll=0; ll<alignmap.size(); ll++) {
626 0 : l=alignmap[ll];
627 0 : int ind=i*3*3*natoms+j*3*natoms+k*natoms+l;
628 0 : tmp1=sqrt(align[l]/totalign);
629 0 : dd_dr_temp[l]=tmp1*dmatdp0[ind];
630 0 : if(do_center) {
631 0 : for(nn=0; nn<alignmap.size(); nn++) {
632 0 : n=alignmap[nn];
633 0 : dd_dr_temp[l]-=dmatdp0[ind-l+n]*tmp1*align[n]/totalign;
634 : }
635 : }
636 : }
637 0 : for(ll=0; ll<alignmap.size(); ll++) {
638 0 : l=alignmap[ll];
639 0 : int ind=i*3*3*natoms+j*3*natoms+k*natoms+l;
640 0 : dmatdp0[ind]=dd_dr_temp[l];
641 : }
642 : }
643 : }
644 : }
645 : }
646 :
647 :
648 0 : bool do_p1rotated=true;
649 0 : if (do_p1rotated) {
650 : // resize if not allocated
651 :
652 0 : if(p1.size()!=p1rotated.size())p1rotated.resize(p1.size());
653 :
654 : // exit(0);
655 :
656 0 : for(i=0; i<natoms; i++) p1rotated[i]=matmul(d,p1reset[i]);
657 :
658 : // reallocate difference vectors
659 0 : if(p1.size()!=diff1on0.size())diff1on0.resize(p1.size());
660 0 : for(i=0; i<natoms; i++) {
661 0 : diff1on0[i]=p1rotated[i]-p0reset[i];
662 : }
663 :
664 0 : if(verbose) {
665 0 : log->printf("P1-RESET-AND-ROTATED\n");
666 0 : for(i=0; i<natoms; i++) {
667 0 : log->printf("ATOM %6u 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 0 : log->printf("END\n");
670 0 : log->printf("P0-RESET\n");
671 0 : for(i=0; i<natoms; i++) {
672 0 : log->printf("ATOM %6u 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 0 : log->printf("END\n");
675 : }
676 :
677 : }
678 :
679 :
680 :
681 0 : dinv=inverse(d);
682 :
683 0 : bool do_p0rotated=true;
684 0 : if (do_p0rotated) {
685 0 : if(p0.size()!=p0rotated.size())p0rotated.resize(p0.size());
686 0 : for(i=0; i<natoms; i++) p0rotated[i]=matmul(dinv,p0reset[i]);
687 0 : if(p1.size()!=diff0on1.size())diff0on1.resize(p1.size());
688 0 : for(i=0; i<natoms; i++) diff0on1[i]=p0rotated[i]-p1reset[i];
689 0 : if(verbose) {
690 0 : log->printf("P0-RESET AND INVERSE ROTATED\n");
691 0 : for(i=0; i<natoms; i++) {
692 0 : log->printf("ATOM %6u 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 0 : log->printf("END\n");
695 0 : log->printf("P1-RESET\n");
696 0 : for(i=0; i<natoms; i++) {
697 0 : log->printf("ATOM %6u 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 0 : log->printf("END\n");
700 : }
701 : }
702 : // copy on the official vectors:
703 0 : rotmat0on1=d;
704 0 : rotmat1on0=dinv;
705 0 : derrdp0.resize(natoms);
706 0 : derrdp0=derr_dr0;
707 0 : derrdp1.resize(natoms);
708 0 : derrdp1=derr_dr1;
709 :
710 : // now rescale accordingly for rmsd instead of msd
711 0 : if(rmsd) {
712 0 : err=sqrt(err);
713 0 : double tmp=0.5/err;
714 0 : for(ii=0; ii<alignmap.size(); ii++) {
715 0 : i=alignmap[ii];
716 0 : derrdp0[i]*=tmp;
717 0 : derrdp1[i]*=tmp;
718 : }
719 :
720 : }
721 :
722 0 : return err;
723 :
724 : }
725 0 : void Kearsley::assignP1(const std::vector<Vector> & p1) {
726 0 : this->p1=p1;
727 0 : com1_is_removed=false;
728 0 : }
729 0 : void Kearsley::assignP0(const std::vector<Vector> & p0) {
730 0 : this->p0=p0;
731 0 : com0_is_removed=false;
732 0 : }
733 :
734 0 : void Kearsley::assignAlign(const std::vector<double> & align) {
735 0 : this->align=align;
736 0 : }
737 :
738 0 : void Kearsley::finiteDifferenceInterface(bool rmsd) {
739 0 : log->printf("Entering rmsd finite difference test system for kearsley\n");
740 0 : log->printf("-------------------------------------------\n");
741 0 : log->printf("TEST1: derivative of the value (derr_dr0/derr_dr1)\n");
742 : //// test 1
743 : unsigned i,j,l,m;
744 0 : double step=1.e-6,olderr,delta;
745 : // messing up a bit with align weights
746 : double delta1;
747 0 : vector<double> align1;
748 0 : align1.resize(p0.size());
749 0 : Random rnd;
750 0 : for (i=0; i<p0.size(); i++) {
751 : // draw a random number
752 0 : delta=rnd.RandU01();
753 0 : delta1=rnd.RandU01();
754 0 : if(delta>delta1) {
755 : //if(delta>0.3){
756 0 : align1[i]=delta;
757 0 : } else {align1[i]=0.;};
758 : // log->printf("ALIGN %d IS %8.3f\n",i,align1[i]);
759 : }
760 0 : assignAlign(align1);
761 : //// get initial value of the error and derivative of it
762 0 : olderr=calculate(rmsd);
763 0 : log->printf("INITIAL ERROR VALUE: %e\n",olderr);
764 : // store the matrix
765 0 : Tensor old_rotmat0on1=rotmat0on1;
766 :
767 : //// get initial value of the error and derivative of it
768 :
769 0 : log->printf("TESTING: derrdp1 \n");
770 0 : for(unsigned j=0; j<3; j++) {
771 0 : for(unsigned i=0; i<derrdp1.size(); i++) {
772 : // random displacement
773 0 : delta=(rnd.RandU01()-0.5)*2*step;
774 0 : p1[i][j]+=delta;
775 0 : com1_is_removed=false; // this is required whenever the assignment is not done with the methods
776 0 : com0_is_removed=false; // this is required whenever the assignment is not done with the methods
777 0 : 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 0 : p1[i][j]-=delta;
780 0 : switch(j) {
781 : case 0:
782 0 : 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 0 : 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 0 : 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 0 : log->printf("TESTING: derrdp0 \n");
793 0 : for(unsigned j=0; j<3; j++) {
794 0 : for(unsigned i=0; i<derrdp0.size(); i++) {
795 : // random displacement
796 0 : delta=(rnd.RandU01()-0.5)*2*step;
797 0 : p0[i][j]+=delta;
798 0 : com0_is_removed=false; // this is required whenever the assignment is not done with the methods
799 0 : com1_is_removed=false; // this is required whenever the assignment is not done with the methods
800 :
801 0 : err=calculate(rmsd);
802 0 : p0[i][j]-=delta;
803 0 : switch(j) {
804 : case 0:
805 0 : 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 0 : 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 0 : 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 0 : log->printf("TESTING: dmatdp0 \n");
815 0 : for(l=0; l<3; l++) {
816 0 : for(m=0; m<3; m++) {
817 0 : for(j=0; j<3; j++) {
818 0 : for(i=0; i<p0.size(); i++) {
819 : // random displacement
820 0 : delta=(rnd.RandU01()-0.5)*2*step;
821 0 : p0[i][j]+=delta;
822 0 : com0_is_removed=false;
823 0 : com1_is_removed=false;
824 0 : calculate(rmsd);
825 0 : p0[i][j]-=delta;
826 0 : int ind=l*3*3*p0.size()+m*3*p0.size()+j*p0.size()+i;
827 0 : switch(j) {
828 : case 0:
829 0 : 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 0 : 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 0 : 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 0 : log->printf("TESTING: dmatdp1 \n");
843 0 : for(l=0; l<3; l++) {
844 0 : for(m=0; m<3; m++) {
845 0 : for(j=0; j<3; j++) {
846 0 : for(i=0; i<p1.size(); i++) {
847 : // random displacement
848 0 : delta=(rnd.RandU01()-0.5)*2*step;
849 0 : p1[i][j]+=delta;
850 0 : com0_is_removed=false;
851 0 : com1_is_removed=false;
852 0 : calculate(rmsd);
853 0 : p1[i][j]-=delta;
854 0 : int ind=l*3*3*p1.size()+m*3*p1.size()+j*p1.size()+i;
855 0 : switch(j) {
856 :
857 : case 0:
858 0 : 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 0 : 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 0 : 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 0 : exit(0);
873 : }
874 2523 : }
|
