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 myR0, double myD, double myBeta) |
23 |
: BondType(myR0), De(myD), beta(myBeta) { |
24 |
} |
25 |
|
26 |
void setWellDepth(double myD) { De = myD;} |
27 |
|
28 |
void setBeta(double myBeta) { beta = myBeta; } |
29 |
|
30 |
void setWellDepthAndForceConstant(double myD, double myK) { |
31 |
De = myD; |
32 |
beta = sqrt(myK/(2.0*De)); |
33 |
} |
34 |
|
35 |
double getWellDepth() {return De;} |
36 |
|
37 |
double getBeta() {return beta;} |
38 |
|
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 |
|
57 |
}; |
58 |
} |
59 |
#endif |