| 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 |
virtual 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 |
