1 |
#ifndef TYPES_MORSEBONDTYPE_HPP |
2 |
#define TYPES_MORSEBONDTYPE_HPP |
3 |
|
4 |
#include "types/BondType.hpp" |
5 |
|
6 |
namespace oopse { |
7 |
/** |
8 |
* @class MorseBondType |
9 |
* |
10 |
* MorseBondType is a more realistic bond potential. |
11 |
* |
12 |
* The functional form is given by: \f$V(r) = D_e (1 - e^{-\beta (r |
13 |
* - r_0)})^2\f$ where \f$D_e\f$ is the bond dissociation energy (in |
14 |
* kcal / mol), \f$\beta$\f is an inverse distance parameter related |
15 |
* to the force constant. \f$\beta = \sqrt{\frac{k}{2 D_e}}\f$, and |
16 |
* \f$r_0\f$ is the equilibrium bond length. |
17 |
*/ |
18 |
class MorseBondType : public BondType { |
19 |
|
20 |
public: |
21 |
|
22 |
MorseBondType(double myD, double myBeta, double myR0) : BondType() { |
23 |
De = myD; |
24 |
beta = myBeta; |
25 |
r0 = myR0; |
26 |
} |
27 |
|
28 |
void setWellDepth(double myD) { De = myD;} |
29 |
void setEquilibriumBondLength(double r) { r0 = r; } |
30 |
void setBeta(double myBeta) { beta = myBeta; } |
31 |
void setWellDepthAndForceConstant(double myD, double myK) { |
32 |
De = myD; |
33 |
beta = sqrt(myK/(2.0*De)); |
34 |
} |
35 |
|
36 |
double getWellDepth() {return De;} |
37 |
double getEquilibriumBondLength() {return r0;} |
38 |
double getBeta() {return beta;} |
39 |
double getForceConstant() {return 2.0*De*beta*beta;} |
40 |
|
41 |
void calcForce(double r, double& V, double& dVdr) { |
42 |
double dr, eterm, eterm2; |
43 |
|
44 |
dr = r - r0; |
45 |
eterm = exp(-beta*dr); |
46 |
eterm2 = eterm*eterm; |
47 |
|
48 |
V = De*(1 - 2.0*eterm + eterm2); |
49 |
dVdr = 2.0 * De * beta * (eterm - eterm2); |
50 |
} |
51 |
|
52 |
private: |
53 |
|
54 |
double De; |
55 |
double beta; |
56 |
double r0; |
57 |
|
58 |
}; |
59 |
} |
60 |
#endif |