| 38 |
|
* University of Notre Dame has been advised of the possibility of |
| 39 |
|
* such damages. |
| 40 |
|
*/ |
| 41 |
+ |
#include <algorithm> |
| 42 |
|
#include <cmath> |
| 43 |
|
#include "types/CharmmTorsionType.hpp" |
| 44 |
|
#include "utils/NumericConstant.hpp" |
| 45 |
+ |
#include "math/ChebyshevPolynomials.hpp" |
| 46 |
|
namespace oopse { |
| 47 |
< |
|
| 46 |
< |
void CharmmTorsionType::calcForce(double cosPhi, double sinPhi, double& V, double& dVdPhi) { |
| 47 |
< |
V = 0; |
| 48 |
< |
dVdPhi = 0; |
| 49 |
< |
double phi = -std::atan2(sinPhi, cosPhi); |
| 50 |
< |
|
| 47 |
> |
CharmmTorsionType::CharmmTorsionType(std::vector<CharmmTorsionParameter>& parameters) { |
| 48 |
|
std::vector<CharmmTorsionParameter>::iterator i; |
| 49 |
< |
for (i = parameter_.begin(); i != parameter_.end(); ++i) { |
| 50 |
< |
double kchi= i->kchi; |
| 51 |
< |
int n = i->n; |
| 52 |
< |
double delta = i->delta; |
| 49 |
> |
i = std::max_element(parameters.begin(), parameters.end(), LessThanPeriodicityFunctor()); |
| 50 |
> |
if (i != parameters.end()) { |
| 51 |
> |
int maxPower = i->n; |
| 52 |
> |
ChebyshevT T(maxPower); |
| 53 |
> |
ChebyshevU U(maxPower); |
| 54 |
|
|
| 55 |
< |
if (n == 0) { |
| 56 |
< |
//if periodicity is equal to 0, use harmonic form |
| 57 |
< |
double diff = phi - delta; |
| 55 |
> |
// convert parameters of charmm type torsion into PolynomialTorsion's parameters |
| 56 |
> |
DoublePolynomial finalPolynomial; |
| 57 |
> |
for (i = parameters.begin(); i != parameters.end(); ++i) { |
| 58 |
> |
DoublePolynomial cosTerm = T.getChebyshevPolynomial(i->n); |
| 59 |
> |
cosTerm *= cos(i->delta) * i->kchi; |
| 60 |
> |
DoublePolynomial sinTerm = U.getChebyshevPolynomial(i->n); |
| 61 |
> |
sinTerm *= -sin(i->delta) * i->kchi; |
| 62 |
> |
finalPolynomial = cosTerm + sinTerm; |
| 63 |
> |
finalPolynomial += i->kchi; |
| 64 |
> |
} |
| 65 |
> |
torsionType_->setPolynomial(finalPolynomial); |
| 66 |
> |
} |
| 67 |
> |
} |
| 68 |
|
|
| 61 |
– |
if (diff < -NumericConstant::PI) { |
| 62 |
– |
diff += NumericConstant::TWO_PI; |
| 63 |
– |
} else if (diff > NumericConstant::PI) { |
| 64 |
– |
diff -= NumericConstant::TWO_PI; |
| 65 |
– |
} |
| 66 |
– |
|
| 67 |
– |
V += kchi * diff * diff; |
| 68 |
– |
dVdPhi += 2.0 * kchi * diff; |
| 69 |
– |
|
| 70 |
– |
} else { |
| 71 |
– |
//use normal cos form if periodicity is greater than 0 |
| 72 |
– |
V += kchi * (1 + std::cos(n * phi + delta)); |
| 73 |
– |
dVdPhi += -n * kchi * std::sin(n * phi + delta); |
| 74 |
– |
} |
| 75 |
– |
} |
| 76 |
– |
} |
| 77 |
– |
|
| 69 |
|
} //end namespace oopse |
