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 "OptimalAlignment.h"
23 : #include "Kearsley.h"
24 : #include "Log.h"
25 : #include <cmath>
26 : #include <iostream>
27 : #include <cstdlib>
28 : #include "Random.h"
29 :
30 : using namespace std;
31 : namespace PLMD {
32 :
33 0 : OptimalAlignment::OptimalAlignment( const std::vector<double> & align, const std::vector<double> & displace, const std::vector<Vector> & p0, const std::vector<Vector> & p1, Log* &log )
34 0 : :log(log) {
35 :
36 : // kearsley init to null
37 0 : mykearsley=NULL;
38 0 : if (mykearsley==NULL) {
39 0 : mykearsley=new Kearsley(p0,p1,align,log);
40 : }
41 : // copy the structure into place
42 0 : assignP0(p0);
43 0 : assignP1(p1);
44 0 : assignAlign(align);
45 0 : assignDisplace(displace);
46 :
47 : // basic check
48 0 : if(p0.size() != p1.size()) {
49 0 : (this->log)->printf("THE SIZE OF THE TWO FRAMES TO BE ALIGNED ARE DIFFERENT\n");
50 : }
51 : // fast behaviour: if all the alignment and displacement are 1.0 then go for fast
52 0 : fast=true;
53 0 : for (unsigned i=0; i<align.size(); i++ ) {
54 0 : if(align[i]!=displace[i] || align[i]!=1.0)fast=false;
55 : }
56 :
57 0 : }
58 :
59 0 : OptimalAlignment::~OptimalAlignment() {
60 0 : if (mykearsley!=NULL)delete mykearsley;
61 0 : }
62 :
63 0 : void OptimalAlignment::assignP0( const std::vector<Vector> & p0 ) {
64 0 : this->p0=p0;
65 0 : if(mykearsley!=NULL) {mykearsley->assignP0(p0);} else {cerr<<"kearsley is not initialized"<<endl; exit(0);}
66 0 : }
67 :
68 0 : void OptimalAlignment::assignP1( const std::vector<Vector> & p1 ) {
69 0 : this->p1=p1;
70 0 : if(mykearsley!=NULL) {mykearsley->assignP1(p1);} else {cerr<<"kearsley is not initialized"<<endl; exit(0);}
71 0 : }
72 :
73 0 : void OptimalAlignment::assignAlign( const std::vector<double> & align ) {
74 0 : this->align=align;
75 0 : if(mykearsley!=NULL) {mykearsley->assignAlign(align);} else {cerr<<"kearsley is not initialized"<<endl; exit(0);}
76 0 : }
77 :
78 0 : void OptimalAlignment::assignDisplace( const std::vector<double> & displace ) {
79 0 : this->displace=displace;
80 0 : }
81 :
82 :
83 0 : double OptimalAlignment::calculate(bool squared, std::vector<Vector> & derivatives) {
84 :
85 0 : bool rmsd=!squared ;
86 :
87 : double err;
88 :
89 : // at this point everything should be already in place for calculating the alignment (p1,p2,align)
90 : // here everything is done with kearsley algorithm. Extension to other optimal alignment algos is
91 : // possible here below with a switch
92 :
93 0 : err=mykearsley->calculate(rmsd); // this calculates the MSD: transform into RMSD
94 :
95 : // check findiff alignment
96 : //mykearsley->finiteDifferenceInterface(rmsd);
97 :
98 0 : if(fast) {
99 : //log->printf("Doing fast: ERR %12.6f \n",err);
100 0 : derrdp0=mykearsley->derrdp0;
101 0 : derrdp1=mykearsley->derrdp1;
102 0 : derivatives=derrdp0;
103 : } else {
104 : /// TODO this interface really sucks since is strongly asymmetric should be re-engineered.
105 0 : err=weightedAlignment(rmsd);
106 : //log->printf("Doing slow: ERR %12.6f \n",err);
107 0 : derivatives=derrdp0;
108 : }
109 : // destroy the kearsley object?
110 :
111 0 : return err;
112 : }
113 :
114 : #ifdef __INTEL_COMPILER
115 : #pragma intel optimization_level 2
116 : #endif
117 : /// this does the weighed alignment if the vector of alignment is different from displacement
118 0 : double OptimalAlignment::weightedAlignment( bool rmsd) {
119 : double tmp0,tmp1,walign,wdisplace,const1,ret;
120 : unsigned i,k,l,m,n,o,oo,mm;
121 :
122 0 : unsigned natoms=p0.size();
123 :
124 0 : Kearsley *myk=mykearsley; /// just a shortcut
125 :
126 : /// TODO : these two blocks can be calculated once forever after the initialization (exception for certain methods?)
127 :
128 : /// clear derivatives
129 0 : if (derrdp0.size()!=natoms)derrdp0.resize(natoms);
130 0 : if (derrdp1.size()!=natoms)derrdp1.resize(natoms);
131 :
132 : // clear the container
133 0 : for(i=0; i<natoms; i++) {
134 0 : derrdp0[i][0]=derrdp0[i][1]=derrdp0[i][2]=0.;
135 0 : derrdp1[i][0]=derrdp1[i][1]=derrdp1[i][2]=0.;
136 : }
137 :
138 0 : walign=0.;
139 0 : vector<int> alignmap;
140 0 : for(i=0; i<natoms; i++) {
141 0 : if (align[i]>0.) {
142 0 : alignmap.push_back(i);
143 0 : walign+=align[i];
144 : }
145 0 : if (align[i]<0.) {cerr<<"FOUND ALIGNMENT WEIGHT NEGATIVE!"<<endl; exit(0);};
146 : }
147 :
148 0 : wdisplace=0.;
149 0 : vector<int> displacemap;
150 0 : for(i=0; i<natoms; i++) {
151 0 : if (displace[i]>0.) {
152 0 : displacemap.push_back(i);
153 0 : wdisplace+=displace[i];
154 : }
155 0 : if (displace[i]<0.) {cerr<<"FOUND ALIGNMENT WEIGHT NEGATIVE!"<<endl; exit(0);};
156 : }
157 :
158 :
159 0 : tmp0=0.;
160 :
161 0 : vector<Vector> array_n_3;
162 0 : array_n_3.resize(natoms);
163 0 : for(i=0; i<array_n_3.size(); i++)array_n_3[i][0]=array_n_3[i][1]=array_n_3[i][2]=0.;
164 :
165 : // err= (1/totdisplace) sum_k_l displace_k*((p0reset_k_l- sum_k_m rot_l_m*p1reset_k_m )**2)
166 :
167 : //for(kk=0;kk<displacemap.size();kk++){
168 : // k=displacemap[kk];
169 0 : for(k=0; k<natoms; k++) {
170 :
171 0 : for(l=0; l<3; l++) {
172 :
173 0 : tmp1=0.;
174 : // contribution from rotated reference frame //
175 0 : for(m=0; m<3; m++) {
176 0 : tmp1-=myk->rotmat0on1[l][m]*myk->p1reset[k][m];
177 : }
178 :
179 : // contribution from running centered frame //
180 0 : tmp1+= myk->p0reset[k][l];
181 :
182 0 : array_n_3[k][l]=tmp1; // store coefficents for derivative usage//
183 0 : tmp0+=tmp1*tmp1*displace[k]; //squared distance added//
184 : }
185 :
186 : }
187 :
188 0 : tmp0=tmp0/wdisplace;
189 :
190 : // log->printf(" ERRR NEW %f \n",tmp0);
191 :
192 0 : ret=tmp0;
193 :
194 : /* DERIVATIVE CALCULATION:respect to running frame */
195 0 : for(k=0; k<natoms; k++) {
196 0 : for(l=0; l<3; l++) {
197 :
198 0 : tmp1 =2.*array_n_3[k][l]*displace[k]/wdisplace ; //ok
199 :
200 0 : const1=2.*align[k]/(walign*wdisplace);
201 :
202 0 : if(const1>0.) {
203 0 : for(oo=0; oo<displacemap.size(); oo++) {
204 0 : o=displacemap[oo];
205 0 : tmp1 -=const1*array_n_3[o][l]*displace[o]; //ok
206 : }
207 : }
208 :
209 0 : for(mm=0; mm<displacemap.size(); mm++) {
210 0 : m=displacemap[mm];
211 0 : const1=2.* displace[m]/wdisplace ;
212 0 : for(n=0; n<3; n++) {
213 0 : tmp0=0.;
214 0 : for(o=0; o<3; o++) {
215 0 : int ind=n*3*3*natoms+o*3*natoms+l*natoms+k; //ok
216 0 : tmp0+=myk->dmatdp0[ind]*myk->p1reset[m][o];
217 : }
218 0 : tmp0*=-const1*array_n_3[m][n];
219 :
220 0 : tmp1+=tmp0;
221 : }
222 : }
223 :
224 0 : derrdp0[k][l]=tmp1;
225 :
226 : }
227 : }
228 : //exit(0);
229 :
230 : //return ret;
231 0 : bool do_frameref_der=true;
232 :
233 : // derivatives of
234 : // err= (1/totdisplace) sum_k_l displace_k*((p0reset_k_l- sum_m rot_l_m*p1reset_k_m )**2)
235 : // respect p1:
236 : // derr_dp1=(1/totdisplace) sum_k_l 2*displace_k*(p0reset_k_l- sum_m rot_l_m*p1reset_k_m )
237 : // *d/dp1 ( p0reset_k_l- sum_m rot_l_m*p1reset_k_m)
238 : // =
239 : // (1/totdisplace) sum_k_l 2*displace_k*(p0reset_k_l- sum_m rot_l_m*p1reset_k_m )*
240 : // *(d/dp1 p0reset_k_l
241 : // - sum_m (d/dp1 rot_l_m)*p1reset_k_m
242 : // - sum_m rot_l_m*(d/dp1 p1reset_k_m ) )
243 : // =
244 : // sum_k_l 2*displace_k/totdisplace* array_n_3_k_l
245 : // *(- sum_m (d/dp1 rot_l_m)*p1reset_k_m
246 : // - sum_m rot_l_m*(d/dp1 p1reset_k_m ) )
247 :
248 0 : if(do_frameref_der) {
249 0 : for(k=0; k<natoms; k++) {
250 : // for(kk=0;kk<displacemap.size();kk++){
251 : // k=displacemap[kk];
252 :
253 0 : for(l=0; l<3; l++) {
254 :
255 0 : tmp1=0.;
256 0 : for(mm=0; mm<displacemap.size(); mm++) {
257 0 : m=displacemap[mm];
258 0 : const1=2.* displace[m]/wdisplace ;
259 0 : for(n=0; n<3; n++) {
260 0 : tmp0=0.;
261 0 : for(o=0; o<3; o++) {
262 0 : int ind=n*3*3*natoms+o*3*natoms+l*natoms+k;
263 0 : tmp0+=myk->dmatdp1[ind]*myk->p1reset[m][o];
264 : }
265 0 : tmp0*=-const1*array_n_3[m][n];
266 0 : tmp1+= tmp0;
267 : }
268 :
269 : }
270 :
271 0 : tmp0=0.;
272 0 : for(o=0; o<3; o++) {
273 0 : tmp0+=array_n_3[k][o]*myk->rotmat0on1[o][l];
274 : }
275 0 : tmp1+=-tmp0*2.*displace[k]/wdisplace;
276 :
277 0 : tmp0=0.;
278 :
279 0 : for(mm=0; mm<displacemap.size(); mm++) {
280 0 : m=displacemap[mm];
281 0 : for(o=0; o<3; o++) {
282 0 : tmp0+=array_n_3[m][o]*myk->rotmat0on1[o][l]*displace[m];
283 : }
284 : }
285 0 : tmp1 += tmp0*2.*align[k]/(walign*wdisplace);
286 :
287 0 : derrdp1[k][l]=tmp1;
288 :
289 : }
290 : }
291 : }
292 :
293 : /// probably not really the way it should be
294 0 : if (rmsd) {
295 0 : ret=sqrt(ret);
296 0 : double tmp=0.5/ret;
297 0 : for(unsigned i=0; i<natoms; i++) {
298 0 : derrdp0[i][0]=derrdp0[i][0]*tmp;
299 0 : derrdp0[i][1]=derrdp0[i][1]*tmp;
300 0 : derrdp0[i][2]=derrdp0[i][2]*tmp;
301 0 : derrdp1[i][0]=derrdp1[i][0]*tmp;
302 0 : derrdp1[i][1]=derrdp1[i][1]*tmp;
303 0 : derrdp1[i][2]=derrdp1[i][2]*tmp;
304 :
305 : }
306 : }
307 :
308 0 : return ret;
309 : }
310 :
311 0 : double OptimalAlignment::weightedFindiffTest( bool rmsd) {
312 :
313 0 : Random rnd;
314 :
315 0 : log->printf("Entering rmsd finite difference test system\n ");
316 0 : log->printf("RMSD OR MSD: %s\n",(rmsd)?"rmsd":"msd");
317 0 : log->printf("-------------------------------------------\n");
318 0 : log->printf("TEST1: derivative of the value (derr_dr0/derr_dr1)\n");
319 : //// test 1
320 0 : double step=1.e-8,olderr,delta,err;
321 0 : vector<Vector> fakederivatives;
322 0 : fakederivatives.resize(p0.size());
323 0 : fast=false;
324 : // randomizing alignments and displacements
325 : /* for (i=0;i<p0.size();i++){
326 : // draw a random number
327 : delta=drand48();
328 : delta1=drand48();
329 : if(delta>delta1){
330 : align[i]=delta;
331 : }else{align[i]=0.;};
332 : delta=drand48();
333 : delta1=drand48();
334 : if(delta>delta1){
335 : displace[i]=delta;
336 : }else{displace[i]=0.;}
337 : }*/
338 : //// get initial value of the error and derivative of it
339 0 : assignAlign(align);
340 0 : assignDisplace(displace);
341 0 : olderr=calculate(rmsd, fakederivatives);
342 :
343 0 : log->printf("INITIAL ERROR VALUE: %e\n",olderr);
344 :
345 : // randomizing alignments and displacements
346 0 : log->printf("TESTING: derrdp0 \n");
347 :
348 0 : for(unsigned j=0; j<3; j++) {
349 0 : for(unsigned i=0; i<derrdp0.size(); i++) {
350 : // random displacement
351 0 : delta=(rnd.RandU01()-0.5)*2*step;
352 0 : p0[i][j]+=delta;
353 0 : assignP0( p0 );
354 0 : err=calculate(rmsd, fakederivatives);
355 0 : p0[i][j]-=delta;
356 0 : assignP0( p0 );
357 0 : switch(j) {
358 : case 0:
359 0 : log->printf("TESTING: X %4u ANAL %18.9f NUMER %18.9f DELTA %18.9f DISP %6.2f ALIGN %6.2f \n",i,derrdp0[i][j],(err-olderr)/delta,derrdp0[i][j]-(err-olderr)/delta,displace[i],align[i]); break;
360 : case 1:
361 0 : log->printf("TESTING: Y %4u ANAL %18.9f NUMER %18.9f DELTA %18.9f DISP %6.2f ALIGN %6.2f \n",i,derrdp0[i][j],(err-olderr)/delta,derrdp0[i][j]-(err-olderr)/delta,displace[i],align[i]); break;
362 : case 2:
363 0 : log->printf("TESTING: Z %4u ANAL %18.9f NUMER %18.9f DELTA %18.9f DISP %6.2f ALIGN %6.2f \n",i,derrdp0[i][j],(err-olderr)/delta,derrdp0[i][j]-(err-olderr)/delta,displace[i],align[i]); break;
364 :
365 : }
366 : }
367 : }
368 :
369 0 : log->printf("TESTING: derrdp1 \n");
370 0 : for(unsigned j=0; j<3; j++) {
371 0 : for(unsigned i=0; i<derrdp1.size(); i++) {
372 : // random displacement
373 0 : delta=(rnd.RandU01()-0.5)*2*step;
374 0 : p1[i][j]+=delta;
375 0 : assignP1( p1 );
376 0 : err=calculate(rmsd, fakederivatives);
377 0 : p1[i][j]-=delta;
378 0 : assignP1( p1 );
379 0 : switch(j) {
380 : case 0:
381 0 : log->printf("TESTING: X %4u ANAL %18.9f NUMER %18.9f DELTA %18.9f DISP %6.2f ALIGN %6.2f \n",i,derrdp1[i][j],(err-olderr)/delta,derrdp1[i][j]-(err-olderr)/delta,displace[i],align[i]); break;
382 : case 1:
383 0 : log->printf("TESTING: Y %4u ANAL %18.9f NUMER %18.9f DELTA %18.9f DISP %6.2f ALIGN %6.2f \n",i,derrdp1[i][j],(err-olderr)/delta,derrdp1[i][j]-(err-olderr)/delta,displace[i],align[i]); break;
384 : case 2:
385 0 : log->printf("TESTING: Z %4u ANAL %18.9f NUMER %18.9f DELTA %18.9f DISP %6.2f ALIGN %6.2f \n",i,derrdp1[i][j],(err-olderr)/delta,derrdp1[i][j]-(err-olderr)/delta,displace[i],align[i]); break;
386 :
387 : }
388 : }
389 : }
390 0 : exit(0);
391 :
392 : // This is to avoid warnings:
393 0 : return 0.0;
394 :
395 : }
396 :
397 :
398 :
399 2523 : }
|
