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