1 // 2 // Gradsolvmat.h 3 // pbsam_xcode 4 // 5 // Created by David Brookes on 6/20/16. 6 // Copyright © 2016 David Brookes. All rights reserved. 7 // 8 9 #ifndef Gradsolvmat_h 10 #define Gradsolvmat_h 11 12 #include <stdio.h> 13 #include "MyMatrix.h" 14 #include "util.h" 15 #include "Solvmat.h" 16 17 /* 18 References: 19 [1] Yap, E., Head-Gordon, T. 2010. JCTC 20 [2] Yap, E., Head-Gordon, T. 2013. JCTC 21 */ 22 23 /* 24 Base class for gradients of ComplexMoleculeMatrix objects 25 */ 26 class GradCmplxMolMat 27 { 28 protected: 29 int wrt_; // with respect to 30 vector<MyMatrix<Ptx> > mat_; // each gradient has three components (use Pt) 31 int p_; 32 int I_; 33 34 35 public: GradCmplxMolMat(int I,int wrt,int ns,int p)36 GradCmplxMolMat(int I, int wrt, int ns, int p) 37 :wrt_(wrt), p_(p), I_(I), mat_(ns, MyMatrix<Ptx> (p, 2*p+1)) 38 { 39 } 40 get_wrt()41 const int get_wrt() const { return wrt_; } get_I()42 const int get_I() const { return I_; } get_p()43 const int get_p() const { return p_; } get_ns()44 const int get_ns() const { return (int) mat_.size(); } get_mat_knm(int k,int n,int m)45 Ptx get_mat_knm(int k, int n, int m) { return mat_[k](n, m+p_); } 46 get_mat_knm_d(int k,int n,int m,int d)47 cmplx get_mat_knm_d(int k, int n, int m, int d) 48 { 49 if ( d == 0 ) return mat_[k](n, m+p_).x(); 50 else if ( d == 1 ) return mat_[k](n, m+p_).y(); 51 return mat_[k](n, m+p_).z(); 52 } set_mat_knm(int k,int n,int m,Ptx val)53 void set_mat_knm(int k, int n, int m, Ptx val) 54 { mat_[k].set_val(n, m+p_, val); } 55 set_mat_knm_d(int k,int n,int m,int d,cmplx val)56 void set_mat_knm_d(int k, int n, int m, int d, cmplx val) 57 { 58 if ( d == 0 ) mat_[k](n, m+p_).set_x(val); 59 else if ( d == 1 ) mat_[k](n, m+p_).set_y(val); 60 else mat_[k](n, m+p_).set_z(val); 61 } 62 get_mat_k(int k)63 MyMatrix<Ptx> get_mat_k(int k) const { return mat_[k]; } set_mat_k(int k,MyMatrix<Ptx> mat)64 void set_mat_k(int k, MyMatrix<Ptx> mat ) 65 { 66 for (int n = 0; n < p_; n++) 67 for (int m = -n; m <= n; m++) 68 mat_[k].set_val(n, m+p_, mat(n, m+p_)); 69 } 70 add_mat_k(int k,MyMatrix<Ptx> mat)71 void add_mat_k(int k, MyMatrix<Ptx> mat ) 72 { 73 for (int n = 0; n < p_; n++) 74 for (int m = -n; m <= n; m++) 75 mat_[k].set_val(n, m+p_, mat_[k](n, m+p_) + mat(n, m+p_)); 76 } 77 78 void reset_mat(); 79 80 friend ostream & operator<<(ostream & fout, GradCmplxMolMat & M) 81 { 82 for (int k = 0; k < M.get_ns(); k++) 83 { 84 fout << "For sphere " << k << endl; 85 // fout << "{"; 86 for (int d = 0; d < 3; d++) 87 { 88 // fout << "{"; 89 fout << " Dim: " << d << endl; 90 for (int n = 0; n < M.get_p(); n++) 91 { 92 for (int m = 0; m <= n; m++) 93 { 94 double real = M.get_mat_knm_d( k, n, m, d).real(); 95 double imag = M.get_mat_knm_d( k, n, m, d).imag(); 96 if(abs(real) < 1e-15 ) real = 0.0; 97 if(abs(imag) < 1e-15 ) imag = 0.0; 98 fout << "(" << setprecision(7)<< real << ", " << imag << ") "; 99 // fout << setprecision(9) <<"{"<< real << ","<<imag<<"},"; 100 } 101 fout << endl; 102 } 103 // fout << "}," ; 104 fout << endl; 105 } 106 // fout << "}," ; 107 fout << endl; 108 } 109 return fout; 110 } 111 print_kmat(int k)112 void print_kmat(int k) 113 { 114 cout << "molecule " << I_ << " For sphere " << k << endl; 115 for (int d = 0; d < 3; d++) 116 { 117 cout << " Dim: " << d << endl; 118 // cout << "{"; 119 for (int n = 0; n < get_p(); n++) 120 { 121 for (int m = 0; m <= n; m++) 122 { 123 double real = get_mat_knm_d( k, n, m, d).real(); 124 double imag = get_mat_knm_d( k, n, m, d).imag(); 125 if(abs(real) < 1e-15 ) real = 0.0; 126 if(abs(imag) < 1e-15 ) imag = 0.0; 127 cout << setprecision(9) << "(" << real << ", " << imag << ") "; 128 // cout<< setprecision(9) <<"{"<< real << ","<<imag<<"},"; 129 } 130 cout << endl; 131 } 132 // cout << "},"; 133 // cout << endl; 134 } 135 cout << endl; 136 } 137 138 }; 139 140 /* 141 Base class for gradients of NumericalMatrix objects 142 */ 143 class GradNumericalMat 144 { 145 protected: 146 int p_; 147 int I_; 148 int wrt_; // with respect to 149 vector<MyMatrix<Ptx> > mat_cmplx_; // each gradient has 3 parts (use Pt) 150 vector<vector<Pt> > mat_; 151 152 153 154 public: GradNumericalMat(int I,int wrt,int ns,int p)155 GradNumericalMat(int I, int wrt, int ns, int p) 156 :p_(p), mat_(ns), mat_cmplx_(ns, MyMatrix<Ptx> (p, 2*p+1)), I_(I), wrt_(wrt) 157 { 158 } 159 get_wrt()160 const int get_wrt() const { return wrt_; } get_I()161 const int get_I() const { return I_; } get_p()162 const int get_p() const { return p_; } get_ns()163 const int get_ns() const { return (int) mat_.size(); } 164 get_mat_kh(int k,int h)165 Pt get_mat_kh(int k, int h) { return mat_[k][h]; } get_mat_knm(int k,int n,int m)166 Ptx get_mat_knm(int k, int n, int m) { return mat_cmplx_[k](n,m+p_); } set_mat_kh(int k,int h,Pt val)167 void set_mat_kh(int k, int h, Pt val) { mat_[k][h] = val; } get_mat_k(int k)168 vector<Pt> get_mat_k(int k) { return mat_[k]; } get_mat()169 vector<vector<Pt> > get_mat() { return mat_; } get_mat_k_len(int k)170 int get_mat_k_len(int k) { return (int)mat_[k].size(); } 171 172 void reset_mat(int k); 173 174 friend ostream & operator<<(ostream & fout, GradNumericalMat & M) 175 { 176 for (int k = 0; k < M.get_ns(); k++) 177 { 178 fout << "For sphere " << k << endl; 179 for (int d = 0; d < 3; d++) 180 { 181 fout << " Dim: " << d << endl; 182 for (int h = 0; h < M.get_mat_k_len(k); h++) 183 { 184 double real; 185 if (d == 0) real = M.get_mat_kh( k, h).x(); 186 else if (d == 1) real = M.get_mat_kh( k, h).y(); 187 else real = M.get_mat_kh( k, h).z(); 188 if(abs(real) < 1e-15 ) real = 0.0; 189 fout << real << ", "; 190 } 191 fout << endl; 192 } 193 } 194 return fout; 195 } 196 print_kmat(int k)197 void print_kmat(int k) 198 { 199 cout << "molecule " << I_ << " For sphere " << k << endl; 200 for (int d = 0; d < 3; d++) 201 { 202 cout << " Dim: " << d << endl; 203 for (int h = 0; h < get_mat_k_len(k); h++) 204 { 205 double real; 206 if (d == 0) real = get_mat_kh( k, h).x(); 207 else if (d == 1) real = get_mat_kh( k, h).y(); 208 else real = get_mat_kh( k, h).z(); 209 if(abs(real) < 1e-15 ) real = 0.0; 210 cout << real << ", "; 211 } 212 cout << endl; 213 } 214 cout << endl; 215 } 216 print_analytical(int k)217 void print_analytical(int k) 218 { 219 cout << "molecule " << I_ << " For sphere " << k << endl; 220 for (int d = 0; d < 3; d++) 221 { 222 cout << " Dim: " << d << endl; 223 for (int n = 0; n < get_p(); n++) 224 { 225 for (int m = 0; m <= n; m++) 226 { 227 double real, imag; 228 if (d == 0) 229 { 230 real = (get_mat_knm( k, n, m).x()).real(); 231 imag = (get_mat_knm( k, n, m).x()).imag(); 232 } else if (d == 1) 233 { 234 real = (get_mat_knm( k, n, m).y()).real(); 235 imag = (get_mat_knm( k, n, m).y()).imag(); 236 } else 237 { 238 real = (get_mat_knm( k, n, m).z()).real(); 239 imag = (get_mat_knm( k, n, m).z()).imag(); 240 } 241 242 if(abs(real) < 1e-15 ) real = 0.0; 243 if(abs(imag) < 1e-15 ) imag = 0.0; 244 cout << setprecision(9)<< "(" << real << ", " << imag << ") "; 245 } 246 cout << endl; 247 } 248 } 249 } 250 251 }; 252 253 // gradient matrices: 254 class GradFMatrix; 255 class GradHMatrix; 256 class GradWHMatrix; 257 class GradWFMatrix; 258 class GradLHNMatrix; 259 class GradLHMatrix; 260 class GradLFMatrix; 261 262 /* 263 Eq. 12a [2] 264 */ 265 class GradWFMatrix : public GradCmplxMolMat 266 { 267 protected: 268 double eps_; 269 double kappa_; 270 271 272 public: 273 GradWFMatrix(int I, int wrt, int ns, int p, 274 double eps_in, double eps_out, 275 double kappa); 276 277 // calculate values for all spheres in this molecule 278 void calc_all_vals(shared_ptr<BaseMolecule> mol, 279 shared_ptr<BesselCalc> bcalc, 280 shared_ptr<GradHMatrix> dH, 281 shared_ptr<GradFMatrix> dF, 282 shared_ptr<GradLHMatrix> dLH, 283 shared_ptr<GradLHNMatrix> dLHN, 284 shared_ptr<GradLFMatrix> dLF); 285 286 // calcualte values for one sphere in this molecule 287 void calc_val_k(int k, shared_ptr<BaseMolecule> mol, 288 vector<double> besseli, 289 vector<double> besselk, 290 shared_ptr<GradHMatrix> dH, 291 shared_ptr<GradFMatrix> dF, 292 shared_ptr<GradLHMatrix> dLH, 293 shared_ptr<GradLHNMatrix> dLHN, 294 shared_ptr<GradLFMatrix> dLF); 295 get_eps()296 const double get_eps() const { return eps_; } 297 }; 298 299 /* 300 Eq. 12b [2] 301 */ 302 class GradWHMatrix : public GradCmplxMolMat 303 { 304 protected: 305 double kappa_; 306 307 public: 308 GradWHMatrix(int I, int wrt, int ns, int p, double kappa); 309 310 void calc_all_vals(shared_ptr<BaseMolecule> mol, 311 shared_ptr<BesselCalc> bcalc, 312 shared_ptr<GradHMatrix> dH, 313 shared_ptr<GradFMatrix> dF, 314 shared_ptr<GradLHMatrix> dLH, 315 shared_ptr<GradLHNMatrix> dLHN, 316 shared_ptr<GradLFMatrix> dLF); 317 318 void calc_val_k(int k, shared_ptr<BaseMolecule> mol, 319 vector<double> besseli, 320 vector<double> besselk, 321 shared_ptr<GradHMatrix> dH, 322 shared_ptr<GradFMatrix> dF, 323 shared_ptr<GradLHMatrix> dLH, 324 shared_ptr<GradLHNMatrix> dLHN, 325 shared_ptr<GradLFMatrix> dLF); 326 }; 327 328 /* 329 Eq. 11a [2] 330 */ 331 class GradFMatrix : public GradCmplxMolMat 332 { 333 public: 334 GradFMatrix(int I, int wrt, int ns, int p); 335 336 void calc_all_vals(shared_ptr<IEMatrix> IE, 337 shared_ptr<GradWFMatrix> dWF); 338 339 void calc_val_k(int k, shared_ptr<IEMatrix> IE, 340 shared_ptr<GradWFMatrix> dWF); 341 342 Ptx calc_df_P(Pt P, int k, shared_ptr<SHCalc> shcalc); 343 344 }; 345 346 /* 347 Eq. 11b [2] 348 */ 349 class GradHMatrix : public GradCmplxMolMat 350 { 351 protected: 352 double kappa_; 353 354 public: 355 GradHMatrix(int I, int wrt, 356 int ns, int p, double kappa); 357 358 void calc_all_vals(shared_ptr<BaseMolecule> mol, 359 shared_ptr<BesselCalc> bcalc, 360 shared_ptr<IEMatrix> IE, 361 shared_ptr<GradWHMatrix> dWH); 362 363 void calc_val_k(int k, 364 vector<double> besseli, 365 shared_ptr<IEMatrix> IE, 366 shared_ptr<GradWHMatrix> dWH); 367 368 // calculate the gradient of h at point P (Eq. S5a) 369 Ptx calc_dh_P(Pt P, int k, vector<double> besseli, 370 shared_ptr<SHCalc> shcalc); 371 get_kappa()372 double get_kappa() const { return kappa_; } 373 374 }; 375 376 /* 377 Eq. 13 [2] 378 */ 379 class GradLFMatrix : public GradNumericalMat 380 { 381 public: 382 GradLFMatrix(int I, int wrt, int ns, int p); 383 384 void init_k(int k, shared_ptr<BaseMolecule> mol, shared_ptr<GradFMatrix> dF, 385 shared_ptr<SHCalc> shcalc, 386 shared_ptr<PreCalcSH> pre_sh, 387 shared_ptr<ExpansionConstants> _expconst, 388 bool no_pre_sh=false); 389 390 void calc_all_vals(shared_ptr<BaseMolecule> mol, vector<int> interpol, 391 shared_ptr<TMatrix> T, shared_ptr<GradFMatrix> dF, 392 shared_ptr<PreCalcSH> pre_sh, 393 bool no_pre_sh=false); 394 395 void calc_val_k(int k, shared_ptr<BaseMolecule> mol, vector<int> interpol, 396 shared_ptr<TMatrix> T, 397 shared_ptr<GradFMatrix> dF, 398 shared_ptr<PreCalcSH> pre_sh, 399 bool no_pre_sh=false); 400 }; 401 402 /* 403 Eq. 13 [2] 404 */ 405 class GradLHMatrix : public GradNumericalMat 406 { 407 protected: 408 double kappa_; 409 410 public: 411 GradLHMatrix(int I, int wrt, int ns, int p, double kappa); 412 413 void init_k(int k, shared_ptr<BaseMolecule> mol, shared_ptr<GradHMatrix> dH, 414 shared_ptr<SHCalc> shcalc, shared_ptr<BesselCalc> bcalc, 415 shared_ptr<PreCalcSH> pre_sh, 416 shared_ptr<ExpansionConstants> _expconst, 417 bool no_pre_sh=false); 418 419 void calc_all_vals(shared_ptr<BaseMolecule> mol, vector<int> interpol, 420 shared_ptr<TMatrix> T, shared_ptr<GradHMatrix> dH, 421 shared_ptr<PreCalcSH> pre_sh, 422 bool no_pre_sh=false); 423 424 void calc_val_k(int k, shared_ptr<BaseMolecule> mol, vector<int> interpol, 425 shared_ptr<TMatrix> T, shared_ptr<GradHMatrix> dH, 426 shared_ptr<PreCalcSH> pre_sh, 427 bool no_pre_sh=false); 428 }; 429 430 /* 431 Eq. 13 [2] 432 */ 433 class GradLHNMatrix : public GradCmplxMolMat 434 { 435 public: 436 GradLHNMatrix(int I, int wrt, int ns, int p); 437 438 void calc_all_vals(shared_ptr<SystemSAM> sys, shared_ptr<TMatrix> T, 439 vector<shared_ptr<GradCmplxMolMat> > gradT_A, 440 vector<shared_ptr<GradHMatrix> > dH); 441 442 void calc_val_k(int k, shared_ptr<SystemSAM> sys, 443 shared_ptr<TMatrix> T, 444 vector<shared_ptr<GradCmplxMolMat> > gradT_A, 445 vector<shared_ptr<GradHMatrix> > dH, 446 double dist_cut = 10.0); 447 448 }; 449 450 #endif /* Gradsolvmat_h */ 451