include/ewalena/numerics/basis_gto.h

Go to the documentation of this file.

// -------------------------------------------------------------------
// @author Toby D. Young
// @version $Id: basis_gto.h 702 2012-02-13 02:08:58Z oneliefleft $
//
// Copyright 2012 namespace ewalena authors. All rights reserved.
//
// Redistribution and use in source and binary forms, with or without
// modification, are permitted provided that the following conditions
// are met:
//
// 1. Redistributions of source code must retain the above copyright
//    notice, this list of conditions and the following disclaimer.
//
// 2. Redistributions in binary form must reproduce the above
//    copyright notice, this list of conditions and the following
//    disclaimer in the documentation and/or other materials provided
//    with the distribution.
//
// THIS SOFTWARE IS PROVIDED BY THE NAMEPSACE EWALENA AUTHORS ``AS
// IS'' AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
// LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS
// FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE
// NAMESPACE EWALENA AUTHORS OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT,
// INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES
// (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR
// SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION)
// HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT,
// STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
// ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED
// OF THE POSSIBILITY OF SUCH DAMAGE.
//
// The views and conclusions contained in the software and
// documentation are those of the authors and should not be
// interpreted as representing official policies, either expressed or
// implied, of the namespace ewalena authors.
// -------------------------------------------------------------------

#include <ewalena/numerics/basis_base.h>

#ifndef __ewalena_basis_gto_h
#define __ewalena_basis_gto_h

#ifdef EWALENA_HAVE_CONFIG_H
#include <ewalena/base/config.h>
#endif

#include <ewalena/base/utilities.h>

namespace ewalena
{
  template <int dim, typename ValueType = double>
    class GaussianTypeOrbital
    : 
    public BasisBase<dim, ValueType> 
    {
    public:
    
    GaussianTypeOrbital () /* const unsigned int n_contractions */
    :
    BasisBase<3, ValueType> (BasisData<ValueType> ("Gaussian-type orbital"))
    {};
    
    ~GaussianTypeOrbital ();
    
    std::string name () const;
    
    void reinit_alpha (const ValueType &alpha);
    
    ValueType value (const ewalena::Point<dim, ValueType> &point) const;
    
    ValueType grad (const ewalena::Point<dim, ValueType> &point) const;
    
    ValueType laplacian (const ewalena::Point<dim, ValueType> &point) const;
    
    ValueType laplacian (const GaussianTypeOrbital<dim, ValueType> &gaussian_i, 
                         const GaussianTypeOrbital<dim, ValueType> &gaussian_j); 
    
    ValueType integral (const GaussianTypeOrbital<dim, ValueType> &gaussian_i, 
                        const GaussianTypeOrbital<dim, ValueType> &gaussian_j); 
    
    ValueType integral (const GaussianTypeOrbital<dim, ValueType> &gaussian_i, 
                        const GaussianTypeOrbital<dim, ValueType> &gaussian_j,
                        const GaussianTypeOrbital<dim, ValueType> &gaussian_k,
                        const GaussianTypeOrbital<dim, ValueType> &gaussian_l); 
    
    ValueType alpha () const;
    
    private:
    
    ValueType __alpha;

    
    }; /* GaussianTypeOrbital */
  
  
  /*-------------- Inline and Other Functions -----------------------*/
  
  template <int dim, typename ValueType>
    inline 
    ValueType
    laplacian (const GaussianTypeOrbital<dim, ValueType> &gaussian_i,
               const GaussianTypeOrbital<dim, ValueType> &gaussian_j) 
  {
    const ValueType ap          = gaussian_i.alpha ()*gaussian_j.alpha ();
    const ValueType as          = gaussian_i.alpha ()+gaussian_j.alpha ();
    const ValueType ratio       = ap/as;
    const ValueType displacment = gaussian_i.origin ().l2_norm () - gaussian_j.origin ().l2_norm ();
    
    ValueType scalar = 0;
    
    switch (dim) 
      { 
      case 3: 
        scalar = 
          ratio * (3-2*ratio*displacment*displacment) * integral (gaussian_i, gaussian_j);
        break; 
        
      default: 
        assert (false); 
      } 
    return scalar;
    
  }
  
  template <int dim, typename ValueType>
    inline 
    ValueType
    integral (const GaussianTypeOrbital<dim, ValueType> &gaussian_i,
              const GaussianTypeOrbital<dim, ValueType> &gaussian_j) 
  {
    const ValueType radius = gaussian_i.origin ().l2_norm (gaussian_j.origin ());
    const ValueType sum    = gaussian_i.alpha () + gaussian_j.alpha (); 
    const ValueType ratio  = (gaussian_i.alpha () * gaussian_j.alpha ()) / sum;
    
    ValueType scalar = 0.;
    
    switch (dim) 
      { 
      case 3: 
        scalar = 
          std::pow (2, 1.5) * std::pow (ratio/sum, 0.75) * 
          std::exp (-ratio*pow (radius, 2)); 
        break; 
        
      default: 
        assert (false); 
      } 
    
    return scalar;
  }
  
  
  template <int dim, typename ValueType>
    inline 
    ValueType
    integral (const GaussianTypeOrbital<dim, ValueType> &gaussian_i,
              const GaussianTypeOrbital<dim, ValueType> &gaussian_j,
              const GaussianTypeOrbital<dim, ValueType> &gaussian_k,
              const GaussianTypeOrbital<dim, ValueType> &gaussian_l) 
  {
    /* First some intermediate variables. */
    const double as1  = gaussian_i.alpha () + gaussian_k.alpha ();
    const double as2  = gaussian_j.alpha () + gaussian_l.alpha ();
    
    const double as   = 
      gaussian_i.alpha () + gaussian_j.alpha () +
      gaussian_k.alpha () + gaussian_l.alpha ();
    
    const double ap   = 
      gaussian_i.alpha () * gaussian_j.alpha () *
      gaussian_k.alpha () * gaussian_l.alpha ();
    
    const double ratio_01 = gaussian_i.alpha () * gaussian_k.alpha () / as1;
    const double ratio_02 = gaussian_j.alpha () * gaussian_l.alpha () / as2;
    
    /* Displacements and centers of the primitive Gaussians. */
    const double Rp     = ( gaussian_i.alpha ()*gaussian_i.origin ().l2_norm () + 
                            gaussian_k.alpha ()*gaussian_k.origin ().l2_norm () )/as1;
    
    const double Rq     = ( gaussian_j.alpha ()*gaussian_j.origin ().l2_norm () + 
                            gaussian_l.alpha ()*gaussian_l.origin ().l2_norm () )/as2;
    
    const double displacement1 = gaussian_i.origin ().l2_norm ()-gaussian_k.origin ().l2_norm ();
    const double displacement2 = gaussian_j.origin ().l2_norm ()-gaussian_l.origin ().l2_norm ();
    
    /* Helper function F0(t) defined for convenience in computing
       nuclear and two-electron integrals. Define F(0) = 1. */
    const double t = (as1*as2/as) * (Rp-Rq) * (Rp-Rq);
    
    ValueType scalar = 0;
    
    switch (dim) 
      { 
      case 3: 
        scalar = (16/std::sqrt (Value::PI)) 
          * std::pow (ap, 0.75) / (as1 * as2 * std::sqrt (as))
          * exp (-ratio_01 * displacement1 * displacement1 
                 -ratio_02 * displacement2 * displacement2)
          * ewalena::Math::incomplete_gamma (t);
        break; 
        
      default: 
        assert (false); 
      } 
    return scalar;
  }
  
} /* namespace ewalena */

#endif  /* __ewalena_basis_gto_h */