|
bem.h00001 // 00002 // bem.h 00003 // 00004 // Copyright (C) 1996 Limit Point Systems, Inc. 00005 // 00006 // Author: Curtis Janssen <cljanss@limitpt.com> 00007 // Maintainer: LPS 00008 // 00009 // This file is part of the SC Toolkit. 00010 // 00011 // The SC Toolkit is free software; you can redistribute it and/or modify 00012 // it under the terms of the GNU Library General Public License as published by 00013 // the Free Software Foundation; either version 2, or (at your option) 00014 // any later version. 00015 // 00016 // The SC Toolkit is distributed in the hope that it will be useful, 00017 // but WITHOUT ANY WARRANTY; without even the implied warranty of 00018 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 00019 // GNU Library General Public License for more details. 00020 // 00021 // You should have received a copy of the GNU Library General Public License 00022 // along with the SC Toolkit; see the file COPYING.LIB. If not, write to 00023 // the Free Software Foundation, 675 Mass Ave, Cambridge, MA 02139, USA. 00024 // 00025 // The U.S. Government is granted a limited license as per AL 91-7. 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 // This represents a solvent by a polarization charge on a dielectric 00042 // boundary surface. 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 // This holds the area associated with each vertex. It is used 00068 // to convert charges to charge densities and back. 00069 double* vertex_area_; 00070 00071 // Given charges compute surface charge density. 00072 void charges_to_surface_charge_density(double *charges); 00073 00074 // Given surface charge density compute charges. 00075 void surface_charge_density_to_charges(double *charges); 00076 public: 00077 BEMSolvent(const Ref<KeyVal>&); 00078 virtual ~BEMSolvent(); 00079 00080 // This should be called after everything is setup--the 00081 // molecule has the correct the geometry and all of the 00082 // parameters have been adjusted. 00083 void init(); 00084 // This gets rid of the system matrix inverse and, optionally, the surface. 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 // NOTE: call allocation routines after init and free routines before done 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 // Given the efield dotted with the normals at the charge positions this 00109 // will compute a new set of charges. 00110 void compute_charges(double* efield_dot_normals, double* charge); 00111 00112 // Given a set of charges and a total charge, this will normalize 00113 // the integrated charge to the charge that would be expected on 00114 // the surface if the given total charge were enclosed within it. 00115 void normalize_charge(double enclosed_charge, double* charges); 00116 00117 // Given charges and nuclear charges compute their interation energy. 00118 double nuclear_charge_interaction_energy(double *nuclear_charge, 00119 double** charge_positions, 00120 double* charge); 00121 00122 // Given charges compute the interaction energy between the nuclei 00123 // and the point charges. 00124 double nuclear_interaction_energy(double** charge_positions, 00125 double* charge); 00126 00127 // Given charges compute the interaction energy for just the surface. 00128 double self_interaction_energy(double** charge_positions, double *charge); 00129 00130 // Given the charges, return the total polarization charge on the surface. 00131 double polarization_charge(double* charge); 00132 00133 // Return the area (available after compute_charges called). 00134 double area() const { return area_; } 00135 // Return the volume (available after compute_charges called). 00136 double volume() const { return volume_; } 00137 // Return the enclosed charge (available after compute_charges called). 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 // this never needs to be called explicitly, but is here now for debugging 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 // Local Variables: 00159 // mode: c++ 00160 // c-file-style: "CLJ" 00161 // End: Generated at Fri Jan 10 08:14:08 2003 for MPQC 2.1.3 using the documentation package Doxygen 1.2.14. |