//
// bem.h
//
// Copyright (C) 1996 Limit Point Systems, Inc.
//
// Author: Curtis Janssen <cljanss@limitpt.com>
// Maintainer: LPS
//
// This file is part of the SC Toolkit.
//
// The SC Toolkit is free software; you can redistribute it and/or modify
// it under the terms of the GNU Library General Public License as published by
// the Free Software Foundation; either version 2, or (at your option)
// any later version.
//
// The SC Toolkit is distributed in the hope that it will be useful,
// but WITHOUT ANY WARRANTY; without even the implied warranty of
// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
// GNU Library General Public License for more details.
//
// You should have received a copy of the GNU Library General Public License
// along with the SC Toolkit; see the file COPYING.LIB.  If not, write to
// the Free Software Foundation, 675 Mass Ave, Cambridge, MA 02139, USA.
//
// The U.S. Government is granted a limited license as per AL 91-7.
//

#ifndef _chemistry_solvent_bem_h
#define _chemistry_solvent_bem_h

#include <util/class/class.h>
#include <util/state/state.h>
#include <util/keyval/keyval.h>
#include <math/isosurf/volume.h>
#include <math/isosurf/surf.h>
#include <math/scmat/matrix.h>
#include <chemistry/molecule/molecule.h>

namespace sc {

// This represents a solvent by a polarization charge on a dielectric
// boundary surface.
class BEMSolvent: public DescribedClass {
  private:
    int debug_;

    Ref<Molecule> solute_;
    Ref<Molecule> solvent_;
    double solvent_density_;
    double dielectric_constant_;
    Ref<SCMatrixKit> matrixkit_;
    RefSCMatrix system_matrix_i_;
    double f_;
    Ref<MessageGrp> grp_;

    double area_;
    double volume_;
    double computed_enclosed_charge_;
    double edisp_;
    double erep_;

    Ref<TriangulatedImplicitSurface> surf_;

    double** alloc_array(int n, int m);
    void free_array(double**);

    // This holds the area associated with each vertex.  It is used
    // to convert charges to charge densities and back.
    double* vertex_area_;

    // Given charges compute surface charge density.
    void charges_to_surface_charge_density(double *charges);

    // Given surface charge density compute charges.
    void surface_charge_density_to_charges(double *charges);
  public:
    BEMSolvent(const Ref<KeyVal>&);
    virtual ~BEMSolvent();

    // This should be called after everything is setup--the
    // molecule has the correct the geometry and all of the
    // parameters have been adjusted.
    void init();
    // This gets rid of the system matrix inverse and, optionally, the surface.
    void done(int clear_surface = 1);

    int ncharge() { return surf_->nvertex(); }

    Ref<Molecule> solvent() { return solvent_ ;}
    double solvent_density() { return solvent_density_ ;}

    // NOTE: call allocation routines after init and free routines before done
    double** alloc_charge_positions() { return alloc_array(ncharge(), 3); }
    void free_charge_positions(double**a) { free_array(a); }

    double** alloc_normals()  { return alloc_array(ncharge(), 3); }
    void free_normals(double**a) { free_array(a); }

    double* alloc_efield_dot_normals()  { return new double[ncharge()]; }
    void free_efield_dot_normals(double*a) { delete[] a; }

    double* alloc_charges() { return new double[ncharge()]; }
    void free_charges(double*a) { delete[] a; }

    void charge_positions(double**);
    void normals(double**);

    // Given the efield dotted with the normals at the charge positions this
    // will compute a new set of charges.
    void compute_charges(double* efield_dot_normals, double* charge);

    // Given a set of charges and a total charge, this will normalize
    // the integrated charge to the charge that would be expected on
    // the surface if the given total charge were enclosed within it.
    void normalize_charge(double enclosed_charge, double* charges);

    // Given charges and nuclear charges compute their interation energy.
    double nuclear_charge_interaction_energy(double *nuclear_charge,
                                             double** charge_positions,
                                             double* charge);

    // Given charges compute the interaction energy between the nuclei
    // and the point charges.
    double nuclear_interaction_energy(double** charge_positions,
                                      double* charge);

    // Given charges compute the interaction energy for just the surface.
    double self_interaction_energy(double** charge_positions, double *charge);
    
    // Given the charges, return the total polarization charge on the surface.
    double polarization_charge(double* charge);

    // Return the area (available after compute_charges called).
    double area() const { return area_; }
    // Return the volume (available after compute_charges called).
    double volume() const { return volume_; }
    // Return the enclosed charge (available after compute_charges called).
    double computed_enclosed_charge() const {
      return computed_enclosed_charge_;
    }

    double disp() {return edisp_;}
    double rep()  {return erep_;}
    double disprep();

    // this never needs to be called explicitly, but is here now for debugging
    void init_system_matrix();

    Ref<TriangulatedImplicitSurface> surface() const { return surf_; }

    Ref<SCMatrixKit> matrixkit() { return matrixkit_; }
};

}

#endif

// Local Variables:
// mode: c++
// c-file-style: "CLJ"
// End: