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