1 //
2 //  PhysCalcCalc.hpp
3 //  pbsam_xcode
4 //
5 //  Created by David Brookes on 6/22/16.
6 //  Copyright © 2016 David Brookes. All rights reserved.
7 //
8 
9 #ifndef PhysCalc_h
10 #define PhysCalc_h
11 
12 #include <stdio.h>
13 #include <memory>
14 #include "Solver.h"
15 #include "BasePhysCalc.h"
16 #include <map>
17 
18 /*
19  Class for calculating interaction energy of a MoleculeSAM given H^(I,k)
20  and LHN^(I,k) matrices
21  */
22 class EnergyCalcSAM : public BaseEnergyCalc
23 {
24 public:
EnergyCalcSAM(int I)25   EnergyCalcSAM(int I) : BaseEnergyCalc(I) { }
26 
27   double calc_energy(shared_ptr<HMatrix> H, shared_ptr<LHNMatrix> LHN);
28 
29   void calc_all_energy(vector<shared_ptr<HMatrix> > H,
30                        vector<shared_ptr<LHNMatrix> > LHN);
31 
32 };
33 
34 class ForceCalcSAM : public BaseForceCalc
35 {
36 protected:
37   shared_ptr<SHCalc> shcalc_;
38   shared_ptr<BesselCalc> bcalc_;
39 
40   // Force on each CG sphere of each molecule
41   vector<vector<Pt> > forces_;
42 
43   int I_; // number of molecules in the system
44   vector<int> ks_;
45 
46   double eps_s_;
47 
48   /*
49    Calculate translational force on a MoleculeSAM given dI_LHN^(I,k), H^(I,k),
50    LHN^(I,k) and dI_H^(I,k)
51    */
52   void calc_fI(int I, shared_ptr<HMatrix> H, shared_ptr<LHNMatrix> LHN,
53                shared_ptr<GradHMatrix> dH, shared_ptr<GradLHNMatrix> dLHN);
54 
55 public:
56 
57   // number of molecules, Ik of all molecules, solvent dielectric:
ForceCalcSAM(int I,vector<int> ki,double es,shared_ptr<SHCalc> shcalc,shared_ptr<BesselCalc> bcalc)58   ForceCalcSAM(int I, vector<int> ki, double es, shared_ptr<SHCalc> shcalc,
59             shared_ptr<BesselCalc> bcalc)
60   :BaseForceCalc(I), shcalc_(shcalc), bcalc_(bcalc), I_(I), ks_(ki), eps_s_(es),
61   forces_(I)
62   {
63   }
64 
65   void calc_all_f(vector<shared_ptr<HMatrix> > H,
66                   vector<shared_ptr<LHNMatrix> > LHN,
67                   vector<vector<shared_ptr<GradHMatrix> > > dH,
68                   vector<vector<shared_ptr<GradLHNMatrix> > > dLHN);
69 
get_all_fIk(int I)70   vector<Pt> get_all_fIk(int I)  {return forces_[I];}
71 };
72 
73 
74 class TorqueCalcSAM : public BaseTorqueCalc
75 {
76 protected:
77   int I_;
78 
79 public:
TorqueCalcSAM(int I)80   TorqueCalcSAM(int I) : BaseTorqueCalc(I), I_(I)  { }
81 
82   void calc_all_tau(shared_ptr<SystemSAM> sys, shared_ptr<ForceCalcSAM> fcalc);
83   Pt calc_tauI(int i, shared_ptr<BaseMolecule> mol, shared_ptr<ForceCalcSAM> fcalc);
84 
85   Pt cross_prod(Pt a, Pt b);
86 };
87 
88 /*
89  Class for calculating energy force and torque in one place
90  */
91 class PhysCalcSAM : public BasePhysCalc
92 {
93 protected:
94   shared_ptr<Solver> _solv_;
95   shared_ptr<GradSolver> _gradSolv_;
96 
97   shared_ptr<EnergyCalcSAM> _eCalc_;
98   shared_ptr<ForceCalcSAM>  _fCalc_;
99   shared_ptr<TorqueCalcSAM> _torCalc_;
100 
101   shared_ptr<SystemSAM> _sys_;
102 
103 public:
104 
105   // constructor just requires an asolver
106   PhysCalcSAM(shared_ptr<Solver> _solv, shared_ptr<GradSolver> _gradsolv,
107            string outfname, Units unit = INTERNAL);
108 
109   void calc_force();
110   void calc_energy();
111   void calc_torque();
112 
113   void print_all();
114 
get_Tau()115   shared_ptr<vector<Pt> > get_Tau() { return _torCalc_->get_tau(); }
get_F()116   shared_ptr<vector<Pt> > get_F()   { return _fCalc_->get_F(); }
get_omega()117   shared_ptr<vector<double> > get_omega() { return _eCalc_->get_omega(); }
118 
get_taui(int i)119   Pt get_taui(int i)        { return _torCalc_->get_taui(i); }
get_forcei(int i)120   Pt get_forcei(int i)      { return _fCalc_->get_forcei(i); }
get_omegai(int i)121   double get_omegai(int i)  { return _eCalc_->get_ei(i); }
get_taui_conv(int i)122   Pt get_taui_conv(int i)        { return _torCalc_->get_taui(i)*unit_conv_; }
get_forcei_conv(int i)123   Pt get_forcei_conv(int i)      { return _fCalc_->get_forcei(i)*unit_conv_; }
get_omegai_conv(int i)124   double get_omegai_conv(int i)  { return _eCalc_->get_ei(i)*unit_conv_; }
get_moli_pos(int i)125   Pt get_moli_pos(int i)    { return _sys_->get_cogi(i); }
126 
127 };
128 
129 #endif /* PhysCalcCalc_h */
130 
131 
132 
133