1 /* +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
2 Copyright (c) 2016-2021 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
23 /*
24 This vast majority of the source code in this file was writting by
25 Sandro Bottaro with some help from Giovanni Bussi
26 */
27
28 #include "ERMSD.h"
29 #include "PDB.h"
30 #include "Matrix.h"
31 #include "Tensor.h"
32
33 #include "Pbc.h"
34 #include <cmath>
35 #include <iostream>
36
37
38 using namespace std;
39 namespace PLMD {
40
41
clear()42 void ERMSD::clear() {
43 reference_mat.clear();
44 }
45
46 //void ERMSD::calcLcs(const vector<Vector> & positions, vector<Vector> &)
47
setReference(const vector<Vector> & reference,const std::vector<unsigned> & pairs_vec,double mycutoff)48 void ERMSD::setReference(const vector<Vector> & reference, const std::vector<unsigned> & pairs_vec, double mycutoff) {
49
50
51 natoms = reference.size();
52 nresidues = natoms/3;
53 unsigned npairs = pairs_vec.size()/2;
54 pairs.resize(npairs);
55 //for(unsigned i=0;i<2*npairs;++i) {
56 // std::cout << "CCC " << pairs_vec[i] << " ";
57 //}
58 for(unsigned i=0; i<npairs; ++i) {
59
60 pairs[i].first = pairs_vec[2*i];
61 pairs[i].second = pairs_vec[2*i+1];
62 }
63
64 cutoff = mycutoff;
65 std::vector<TensorGeneric<4,3> > deri;
66 deri.resize(natoms*natoms);
67 reference_mat.resize(nresidues*nresidues);
68 Pbc fake_pbc;
69
70
71 calcMat(reference,fake_pbc,reference_mat,deri);
72
73 }
74
75
inPair(unsigned i,unsigned j)76 bool ERMSD::inPair(unsigned i, unsigned j) {
77
78 //return true;
79 if(pairs.size()==0) return true;
80 for(unsigned idx=0; idx<pairs.size(); idx++) {
81 //std::cout << "AAA " << pairs[idx][0] << " " << pairs[idx][1] << "\n";
82 if(pairs[idx].first == i && pairs[idx].second == j) return true;
83 if(pairs[idx].second == i && pairs[idx].first == j) return true;
84 }
85 return false;
86 }
87 //double ERMSD::calculate(const std::vector<Vector> & positions,
88 // std::vector<Vector> &derivatives, Tensor& virial) const {
89
90 // Pbc fake_pbc;
91 // return ERMSD::calculate(positions,fake_pbc,derivatives,virial,false);
92 // }
93
94
95
calcMat(const std::vector<Vector> & positions,const Pbc & pbc,std::vector<Vector4d> & mat,std::vector<TensorGeneric<4,3>> & Gderi)96 void ERMSD::calcMat(const std::vector<Vector> & positions,const Pbc& pbc, std::vector<Vector4d> &mat, std::vector<TensorGeneric<4,3> > &Gderi) {
97
98 std::vector<Vector3d> pos;
99 pos.resize(3*nresidues);
100
101 std::vector<Tensor3d> deri;
102 deri.resize(nresidues*9);
103
104 std::vector<Vector> centers;
105 centers.resize(nresidues);
106
107 unsigned idx_deri = 0;
108
109 Tensor da_dxa = (2./3.)*Tensor::identity();
110 Tensor da_dxb = -(1./3.)*Tensor::identity();
111 Tensor da_dxc = -(1./3.)*Tensor::identity();
112
113 Tensor db_dxa = -(1./3.)*Tensor::identity();
114 Tensor db_dxb = (2./3.)*Tensor::identity();
115 Tensor db_dxc = -(1./3.)*Tensor::identity();
116
117 // Form factors - should this be somewhere else?
118
119 double w = 1./3.;
120 Vector form_factor = Vector(2.0,2.0,1.0/0.3);
121
122 for(unsigned res_idx=0; res_idx<natoms/3; res_idx++) {
123
124
125 const unsigned at_idx = 3*res_idx;
126 //center
127 for (unsigned j=0; j<3; j++) {
128 centers[res_idx] += w*positions[at_idx+j];
129 }
130
131 Vector3d a = delta(centers[res_idx],positions[at_idx]);
132 Vector3d b = delta(centers[res_idx],positions[at_idx+1]);
133 Vector3d d = crossProduct(a,b);
134 double ianorm = 1./a.modulo();
135 double idnorm = 1./d.modulo();
136
137 // X vector: COM-C2
138 pos[at_idx] = a*ianorm;
139 // Z versor: C2 x (COM-C4/C6)
140 pos[at_idx+2] = d*idnorm;
141 // Y versor: Z x Y
142 pos[at_idx+1] = crossProduct(pos[at_idx+2],pos[at_idx]);
143
144 // Derivatives ////////
145 Tensor3d t1 = ianorm*(Tensor::identity()-extProduct(pos[at_idx],pos[at_idx]));
146 // dv1/dxa
147 deri[idx_deri] = (2./3. )*t1;
148 // dv1/dxb
149 deri[idx_deri+3] = -(1./3.)*t1;
150 // dv1/dxc
151 deri[idx_deri+6] = -(1./3.)*t1;
152
153 Tensor dd_dxa = VcrossTensor(a,db_dxa) -VcrossTensor(b,da_dxa);
154 Tensor dd_dxb = VcrossTensor(a,db_dxb)-VcrossTensor(b,da_dxb);
155 Tensor dd_dxc = VcrossTensor(a,db_dxc)-VcrossTensor(b,da_dxc);
156
157 // dv3/dxa
158 deri[idx_deri+2] = deriNorm(d,dd_dxa);
159 // dv3/dxb
160 deri[idx_deri+5] = deriNorm(d,dd_dxb);
161 // dv3/dxc
162 deri[idx_deri+8] = deriNorm(d,dd_dxc);
163
164 // dv2/dxa = dv3/dxa cross v1 + v3 cross dv1/dxa
165 deri[idx_deri+1] = (VcrossTensor(deri[idx_deri+2],pos[at_idx]) + \
166 VcrossTensor(pos[at_idx+2],deri[idx_deri]));
167 // dv2/dxb
168 deri[idx_deri+4] = (VcrossTensor(deri[idx_deri+5],pos[at_idx]) + \
169 VcrossTensor(pos[at_idx+2],deri[idx_deri+3]));
170 // dv2/dxc
171 deri[idx_deri+7] = (VcrossTensor(deri[idx_deri+8],pos[at_idx]) + \
172 VcrossTensor(pos[at_idx+2],deri[idx_deri+6]));
173
174 idx_deri += 9;
175 // End derivatives ///////
176
177 }
178
179
180 // Initialization (unnecessary?)
181 for (unsigned i1=0; i1<nresidues*nresidues; i1++) {
182 for (unsigned i2=0; i2<4; i2++) {
183 mat[i1][i2] = 0.0;
184 }
185 }
186
187 double maxdist = cutoff/form_factor[0];
188 double gamma = pi/cutoff;
189 unsigned idx;
190 unsigned idx1 = 0;
191 // Calculate mat
192 for (unsigned i=0; i<nresidues; i++) {
193 for (unsigned j=0; j<nresidues; j++) {
194
195 // skip i==j
196 if(inPair(i,j) and i != j) {
197 //if(i!=j){
198
199
200 // Calculate normal distance first
201 Vector diff = delta(centers[i],centers[j]);
202 double d1 = diff.modulo();
203 //std::cout << inPair(i,j) << " " << i << " " << j << " "<< d1 <<"\n";
204 //std::cout << inPair(i,j) << " " << i << " " << j << " "<< d1 <<"\n";
205 if(d1<maxdist) {
206
207 // calculate r_tilde_ij
208 Vector3d rtilde;
209 for (unsigned k=0; k<3; k++) {
210 for (unsigned l=0; l<3; l++) {
211 rtilde[l] += pos[3*i+l][k]*diff[k]*form_factor[l];
212 }
213 }
214 double rtilde_norm = rtilde.modulo();
215
216 double irnorm = 1./rtilde_norm;
217
218 // ellipsoidal cutoff
219 if(rtilde_norm < cutoff) {
220 idx = i*nresidues + j;
221 //std::cout << i << " " << j << " " << rtilde_norm << " " << idx <<"\n";
222
223
224 // fill 4d matrix
225 double dummy = sin(gamma*rtilde_norm)/(rtilde_norm*gamma);
226 mat[idx][0] = dummy*rtilde[0];
227 mat[idx][1] = dummy*rtilde[1];
228 mat[idx][2] = dummy*rtilde[2];
229 mat[idx][3] = (1.+ cos(gamma*rtilde_norm))/gamma;
230
231 // Derivative (drtilde_dx)
232 std::vector<Tensor3d> drtilde_dx;
233 drtilde_dx.resize(6);
234 unsigned pos_idx = 3*i;
235 unsigned deri_idx = 9*i;
236 for (unsigned at=0; at<3; at++) {
237 for (unsigned l=0; l<3; l++) {
238 Vector3d rvec = form_factor[l]*((pos[pos_idx+l])/3.);
239 Vector3d vvec = form_factor[l]*(matmul(deri[deri_idx+3*at+l],diff));
240 drtilde_dx[at].setRow(l,vvec-rvec);
241 drtilde_dx[at+3].setRow(l,rvec);
242 }
243 }
244
245 //std::vector<TensorGeneric<4,3> > dG_dx;
246 //dG_dx.resize(6);
247
248 double dummy1 = (cos(gamma*rtilde_norm) - dummy);
249
250 idx1 = i*nresidues*6 + j*6;
251
252 for (unsigned l=0; l<6; l++) {
253 //std::cout << i << " " << j << " " << idx1 << " " << idx1+l << "\n";
254
255 // components 1,2,3
256 // sin(gamma*|rtilde|)/gamma*|rtilde|*d_rtilde +
257 // + ((d_rtilde*r_tilde/r_tilde^2) out r_tilde)*
258 // (cos(gamma*|rtilde| - sin(gamma*|rtilde|)/gamma*|rtilde|))
259 Vector3d rdr = matmul(rtilde,drtilde_dx[l]);
260 Tensor tt = dummy*drtilde_dx[l] + (dummy1*irnorm*irnorm)*Tensor(rtilde,rdr);
261 for (unsigned m=0; m<3; m++) {
262 // Transpose here
263 //dG_dx[l].setRow(m,tt.getRow(m));
264 Gderi[idx1+l].setRow(m,tt.getRow(m));
265 }
266 // component 4
267 // - sin(gamma*|rtilde|)/|rtilde|*(r_tilde*d_rtilde)
268 //dG_dx[l].setRow(3,-dummy*gamma*rdr);
269 Gderi[idx1+l].setRow(3,-dummy*gamma*rdr);
270 }
271
272
273
274
275 }
276 }
277 }
278
279 }
280 }
281
282 }
283
284
calculate(const std::vector<Vector> & positions,const Pbc & pbc,std::vector<Vector> & derivatives,Tensor & virial)285 double ERMSD::calculate(const std::vector<Vector> & positions,const Pbc& pbc,\
286 std::vector<Vector> &derivatives, Tensor& virial) {
287
288
289 double ermsd=0.;
290 std::vector<Vector4d> mat;
291 mat.resize(nresidues*nresidues);
292
293 std::vector<TensorGeneric<4,3> > Gderi;
294 Gderi.resize(natoms*natoms);
295
296 calcMat(positions,pbc,mat,Gderi);
297
298 unsigned idx1 = 0;
299 for(unsigned i=0; i<nresidues; i++) {
300 for(unsigned j=0; j<nresidues; j++) {
301 unsigned ii = i*nresidues + j;
302
303 Vector4d dd = delta(reference_mat[ii],mat[ii]);
304 double val = dd.modulo2();
305 //std::cout << "AAA " << i << " " << j << " " << ii << " "<< val << "\n";
306
307 if(val>0.0 && i!=j) {
308
309 for(unsigned k=0; k<3; k++) {
310 idx1 = i*nresidues*6 + j*6 + k;
311
312 derivatives[3*i+k] += matmul(dd,Gderi[idx1]);
313 derivatives[3*j+k] += matmul(dd,Gderi[idx1+3]);
314 }
315 ermsd += val;
316 }
317 }
318 }
319
320 //std::cout << ermsd << " ";
321 //if(pairs.size()!=0) nresidues=pairs.size();
322 //std::cout << ermsd << " " << nresidues;
323 ermsd = sqrt(ermsd/nresidues);
324 double iermsd = 1.0/(ermsd*nresidues);
325 for(unsigned i=0; i<natoms; ++i) {derivatives[i] *= iermsd;}
326
327 return ermsd;
328 }
329
330
331 }
332