1 |
< |
/* |
1 |
> |
/* |
2 |
|
* Copyright (c) 2005 The University of Notre Dame. All Rights Reserved. |
3 |
|
* |
4 |
|
* The University of Notre Dame grants you ("Licensee") a |
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; |
58 |
< |
|
59 |
< |
if (diff < -NumericConstant::PI) { |
60 |
< |
diff += NumericConstant::TWO_PI; |
61 |
< |
} else if (diff > NumericConstant::PI) { |
62 |
< |
diff -= NumericConstant::TWO_PI; |
63 |
< |
} |
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); |
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 |
< |
} |
65 |
> |
torsionType_->setPolynomial(finalPolynomial); |
66 |
> |
} |
67 |
|
} |
68 |
|
|
69 |
|
} //end namespace oopse |