1 // 2 // bem.h 3 // 4 // Copyright (C) 1996 Limit Point Systems, Inc. 5 // 6 // Author: Curtis Janssen <cljanss@limitpt.com> 7 // Maintainer: LPS 8 // 9 // This file is part of the SC Toolkit. 10 // 11 // The SC Toolkit is free software; you can redistribute it and/or modify 12 // it under the terms of the GNU Library General Public License as published by 13 // the Free Software Foundation; either version 2, or (at your option) 14 // any later version. 15 // 16 // The SC Toolkit is distributed in the hope that it will be useful, 17 // but WITHOUT ANY WARRANTY; without even the implied warranty of 18 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 19 // GNU Library General Public License for more details. 20 // 21 // You should have received a copy of the GNU Library General Public License 22 // along with the SC Toolkit; see the file COPYING.LIB. If not, write to 23 // the Free Software Foundation, 675 Mass Ave, Cambridge, MA 02139, USA. 24 // 25 // The U.S. Government is granted a limited license as per AL 91-7. 26 // 27 28 #ifndef _chemistry_solvent_bem_h 29 #define _chemistry_solvent_bem_h 30 31 #include <util/class/class.h> 32 #include <util/state/state.h> 33 #include <util/keyval/keyval.h> 34 #include <math/isosurf/volume.h> 35 #include <math/isosurf/surf.h> 36 #include <math/scmat/matrix.h> 37 #include <chemistry/molecule/molecule.h> 38 39 namespace sc { 40 41 // This represents a solvent by a polarization charge on a dielectric 42 // boundary surface. 43 class BEMSolvent: public DescribedClass { 44 private: 45 int debug_; 46 47 Ref<Molecule> solute_; 48 Ref<Molecule> solvent_; 49 double solvent_density_; 50 double dielectric_constant_; 51 Ref<SCMatrixKit> matrixkit_; 52 RefSCMatrix system_matrix_i_; 53 double f_; 54 Ref<MessageGrp> grp_; 55 56 double area_; 57 double volume_; 58 double computed_enclosed_charge_; 59 double edisp_; 60 double erep_; 61 62 Ref<TriangulatedImplicitSurface> surf_; 63 64 double** alloc_array(int n, int m); 65 void free_array(double**); 66 67 // This holds the area associated with each vertex. It is used 68 // to convert charges to charge densities and back. 69 double* vertex_area_; 70 71 // Given charges compute surface charge density. 72 void charges_to_surface_charge_density(double *charges); 73 74 // Given surface charge density compute charges. 75 void surface_charge_density_to_charges(double *charges); 76 public: 77 BEMSolvent(const Ref<KeyVal>&); 78 virtual ~BEMSolvent(); 79 80 // This should be called after everything is setup--the 81 // molecule has the correct the geometry and all of the 82 // parameters have been adjusted. 83 void init(); 84 // This gets rid of the system matrix inverse and, optionally, the surface. 85 void done(int clear_surface = 1); 86 ncharge()87 int ncharge() { return surf_->nvertex(); } 88 solvent()89 Ref<Molecule> solvent() { return solvent_ ;} solvent_density()90 double solvent_density() { return solvent_density_ ;} 91 92 // NOTE: call allocation routines after init and free routines before done alloc_charge_positions()93 double** alloc_charge_positions() { return alloc_array(ncharge(), 3); } free_charge_positions(double ** a)94 void free_charge_positions(double**a) { free_array(a); } 95 alloc_normals()96 double** alloc_normals() { return alloc_array(ncharge(), 3); } free_normals(double ** a)97 void free_normals(double**a) { free_array(a); } 98 alloc_efield_dot_normals()99 double* alloc_efield_dot_normals() { return new double[ncharge()]; } free_efield_dot_normals(double * a)100 void free_efield_dot_normals(double*a) { delete[] a; } 101 alloc_charges()102 double* alloc_charges() { return new double[ncharge()]; } free_charges(double * a)103 void free_charges(double*a) { delete[] a; } 104 105 void charge_positions(double**); 106 void normals(double**); 107 108 // Given the efield dotted with the normals at the charge positions this 109 // will compute a new set of charges. 110 void compute_charges(double* efield_dot_normals, double* charge); 111 112 // Given a set of charges and a total charge, this will normalize 113 // the integrated charge to the charge that would be expected on 114 // the surface if the given total charge were enclosed within it. 115 void normalize_charge(double enclosed_charge, double* charges); 116 117 // Given charges and nuclear charges compute their interation energy. 118 double nuclear_charge_interaction_energy(double *nuclear_charge, 119 double** charge_positions, 120 double* charge); 121 122 // Given charges compute the interaction energy between the nuclei 123 // and the point charges. 124 double nuclear_interaction_energy(double** charge_positions, 125 double* charge); 126 127 // Given charges compute the interaction energy for just the surface. 128 double self_interaction_energy(double** charge_positions, double *charge); 129 130 // Given the charges, return the total polarization charge on the surface. 131 double polarization_charge(double* charge); 132 133 // Return the area (available after compute_charges called). area()134 double area() const { return area_; } 135 // Return the volume (available after compute_charges called). volume()136 double volume() const { return volume_; } 137 // Return the enclosed charge (available after compute_charges called). computed_enclosed_charge()138 double computed_enclosed_charge() const { 139 return computed_enclosed_charge_; 140 } 141 disp()142 double disp() {return edisp_;} rep()143 double rep() {return erep_;} 144 double disprep(); 145 146 // this never needs to be called explicitly, but is here now for debugging 147 void init_system_matrix(); 148 surface()149 Ref<TriangulatedImplicitSurface> surface() const { return surf_; } 150 matrixkit()151 Ref<SCMatrixKit> matrixkit() { return matrixkit_; } 152 }; 153 154 } 155 156 #endif 157 158 // Local Variables: 159 // mode: c++ 160 // c-file-style: "CLJ" 161 // End: 162