# bem.h

    00001 //
    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: