00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024
00025
00026
00027
00028
00029
00038 #ifndef MOLECULE_HEADER
00039 #define MOLECULE_HEADER
00040
00041 #include <cmath>
00042 #include <vector>
00043 #include <cassert>
00044
00045 #include "realtype.h"
00046 #include "mat_gblas.h"
00047
00051 struct Atom {
00052 ergo_real charge;
00053 ergo_real coords[3];
00054 };
00055
00060 struct Vector3D {
00061 ergo_real v[3];
00062 Vector3D() {}
00063 Vector3D(ergo_real x, ergo_real y, ergo_real z) {
00064 v[0] = x; v[1] = y; v[2] = z;
00065 }
00066 ergo_real& operator[](unsigned i) { return v[i]; }
00067 ergo_real operator[](unsigned i) const { return v[i]; }
00069 ergo_real dist2(const ergo_real b[]) const {
00070 ergo_real d, r;
00071 d = v[0]-b[0]; r = d*d;
00072 d = v[1]-b[1]; r += d*d;
00073 d = v[2]-b[2]; r += d*d;
00074 return r;
00075 }
00077 ergo_real dist(const Vector3D& b) const
00078 { return template_blas_sqrt(dist2(b.v)); }
00079 ergo_real dist(const ergo_real b[]) const
00080 { return template_blas_sqrt(dist2(b)); }
00081 };
00082
00087 class Molecule {
00088 private:
00089 std::vector<Atom> atoms;
00090 ergo_real netCharge;
00091 int noOfAtoms;
00092
00093 public:
00094
00095 Molecule() : atoms(10), netCharge(0), noOfAtoms(0) {}
00096
00097 void addAtom(ergo_real c, ergo_real x, ergo_real y, ergo_real z) {
00098 int currListSize = atoms.size();
00099 if(noOfAtoms >= currListSize)
00100 atoms.resize(currListSize*2);
00101 atoms[noOfAtoms].charge = c;
00102 atoms[noOfAtoms].coords[0] = x;
00103 atoms[noOfAtoms].coords[1] = y;
00104 atoms[noOfAtoms].coords[2] = z;
00105 noOfAtoms++;
00106 }
00107
00108 void clear() { noOfAtoms = 0; netCharge = 0; }
00109 void setNetCharge(ergo_real netCharge_) { netCharge = netCharge_; }
00110 void replaceAtom(int i, const Atom & atom) { assert(i >= 0 && i < noOfAtoms); atoms[i] = atom; }
00111 void setAtomList(const std::vector<Atom> atomList) { atoms = atomList; noOfAtoms = atomList.size(); }
00112 const Atom* getAtomListPtr() const { return &atoms[0]; }
00113 const Atom & getAtom(int i) const { return atoms[i]; }
00114 int getNoOfAtoms() const { return noOfAtoms; }
00115 ergo_real getNetCharge() const { return netCharge; }
00116
00118 void getExtremeInternuclearDistancesQuadratic(ergo_real & minDist, ergo_real & maxDist) const;
00120 ergo_real getNuclearRepulsionEnergyQuadratic() const;
00122 ergo_real getNuclearElectricFieldEnergy(const Vector3D& electricField) const;
00125 int getNumberOfElectrons() const;
00127 void getNuclearRepulsionEnergyGradientContribQuadratic(ergo_real* resultGradient) const;
00128
00132 int setFromMoleculeFile(const char* fileName,
00133 int netCharge,
00134 char **basissetFile);
00135
00136 };
00137
00138
00139 #endif