spce_ccl.hpp

Go to the documentation of this file.

/*
** This file is part of eOn.
**
** SPDX-License-Identifier: BSD-3-Clause
**
** Copyright (c) 2010--present, eOn Development Team
** All rights reserved.
**
** Repo:
** https://github.com/TheochemUI/eOn
*/
#pragma once
#include "ccl.hpp"
 
namespace forcefields {
class SpceCcl : public Ccl {
public:
  SpceCcl();
  SpceCcl(double cutoff, double switchingWidth);

  void computeHH_O_(const int nAtoms, const double R[], double F[], double &U,
                    const double b[]);

  void computeHH_O_(const int nAtoms, const double R[], double F[], double &U,
                    const double b[], const bool fixed[]);
 
  char const *getName() const; 
 
protected:
  struct Water {
    Water(double const rh1[], double const rh2[], double const ro[],
          double const rc[], double fh1[], double fh2[], double fo[])
        : rh1_(rh1), rh2_(rh2), ro_(ro), rc_(rc), fh1_(fh1), fh2_(fh2),
          fo_(fo) {}
    Water(Water const &w, double fh1[], double fh2[], double fo[])
        : rh1_(w.rh1_), rh2_(w.rh2_), ro_(w.ro_), rc_(w.rc_), fh1_(fh1),
          fh2_(fh2), fo_(fo) {}
    void addForces(Water const &w) {
      for (int a = 0; a < 3; ++a) {
        fh1_[a] += w.fh1_[a];
        fh2_[a] += w.fh2_[a];
        fo_[a] += w.fo_[a];
      };
    }
    const double *const rh1_;
    const double *const rh2_;
    const double *const ro_;
    const double *const rc_;
    double *const fh1_;
    double *const fh2_;
    double *const fo_;
  };

  void intramolecular(Water &water, double &U);

  void lennardJonesWithCutoff(Water &w1, Water &w2, double &U);

  void coulombWithCutoff(Water &w1, Water &w2, double &U);

  void coulombFull(Water &w1, Water &w2, double &U);

  static const double roh_;
  static const double theta_;
  static const double rhh_;
  static const double charge_;
  static const double charge2_;
  static const double A_;
  static const double B_;
  static const double sigma_;
  static const double epsilon_;
  static double const polarisationEnergy_;
 
private:
  void initialiseRho(Vector3 const &v, Rho &r);

  template <int H, int O>
  void computeTemplate(const int nMolecules, const double (*const rh1)[H * 3],
                       const double (*const rh2)[H * 3],
                       const double (*const ro)[O * 3],
                       double (*const fh1)[H * 3], double (*const fh2)[H * 3],
                       double (*const fo)[O * 3], double &energy,
                       double const b[], bool const (*const xh1)[H] = 0,
                       bool const (*const xh2)[H] = 0,
                       bool const (*const xo)[O] = 0);
};
} // namespace forcefields