Google

Main Page   Class Hierarchy   Compound List   File List   Compound Members   Related Pages  

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:

Generated at Fri Jan 10 08:14:08 2003 for MPQC 2.1.3 using the documentation package Doxygen 1.2.14.