00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024
00025
00026
00027
00028 #ifndef _chemistry_solvent_bem_h
00029 #define _chemistry_solvent_bem_h
00030
00031 #include <util/class/class.h>
00032 #include <util/state/state.h>
00033 #include <util/keyval/keyval.h>
00034 #include <math/isosurf/volume.h>
00035 #include <math/isosurf/surf.h>
00036 #include <math/scmat/matrix.h>
00037 #include <chemistry/molecule/molecule.h>
00038
00039 namespace sc {
00040
00041
00042
00043 class BEMSolvent: public DescribedClass {
00044 private:
00045 int debug_;
00046
00047 Ref<Molecule> solute_;
00048 Ref<Molecule> solvent_;
00049 double solvent_density_;
00050 double dielectric_constant_;
00051 Ref<SCMatrixKit> matrixkit_;
00052 RefSCMatrix system_matrix_i_;
00053 double f_;
00054 Ref<MessageGrp> grp_;
00055
00056 double area_;
00057 double volume_;
00058 double computed_enclosed_charge_;
00059 double edisp_;
00060 double erep_;
00061
00062 Ref<TriangulatedImplicitSurface> surf_;
00063
00064 double** alloc_array(int n, int m);
00065 void free_array(double**);
00066
00067
00068
00069 double* vertex_area_;
00070
00071
00072 void charges_to_surface_charge_density(double *charges);
00073
00074
00075 void surface_charge_density_to_charges(double *charges);
00076 public:
00077 BEMSolvent(const Ref<KeyVal>&);
00078 virtual ~BEMSolvent();
00079
00080
00081
00082
00083 void init();
00084
00085 void done(int clear_surface = 1);
00086
00087 int ncharge() { return surf_->nvertex(); }
00088
00089 Ref<Molecule> solvent() { return solvent_ ;}
00090 double solvent_density() { return solvent_density_ ;}
00091
00092
00093 double** alloc_charge_positions() { return alloc_array(ncharge(), 3); }
00094 void free_charge_positions(double**a) { free_array(a); }
00095
00096 double** alloc_normals() { return alloc_array(ncharge(), 3); }
00097 void free_normals(double**a) { free_array(a); }
00098
00099 double* alloc_efield_dot_normals() { return new double[ncharge()]; }
00100 void free_efield_dot_normals(double*a) { delete[] a; }
00101
00102 double* alloc_charges() { return new double[ncharge()]; }
00103 void free_charges(double*a) { delete[] a; }
00104
00105 void charge_positions(double**);
00106 void normals(double**);
00107
00108
00109
00110 void compute_charges(double* efield_dot_normals, double* charge);
00111
00112
00113
00114
00115 void normalize_charge(double enclosed_charge, double* charges);
00116
00117
00118 double nuclear_charge_interaction_energy(double *nuclear_charge,
00119 double** charge_positions,
00120 double* charge);
00121
00122
00123
00124 double nuclear_interaction_energy(double** charge_positions,
00125 double* charge);
00126
00127
00128 double self_interaction_energy(double** charge_positions, double *charge);
00129
00130
00131 double polarization_charge(double* charge);
00132
00133
00134 double area() const { return area_; }
00135
00136 double volume() const { return volume_; }
00137
00138 double computed_enclosed_charge() const {
00139 return computed_enclosed_charge_;
00140 }
00141
00142 double disp() {return edisp_;}
00143 double rep() {return erep_;}
00144 double disprep();
00145
00146
00147 void init_system_matrix();
00148
00149 Ref<TriangulatedImplicitSurface> surface() const { return surf_; }
00150
00151 Ref<SCMatrixKit> matrixkit() { return matrixkit_; }
00152 };
00153
00154 }
00155
00156 #endif
00157
00158
00159
00160
00161