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