1 // $Id$
2 //
3 // Created by Santosh Putta, Nov 2006
4 //
5
6 #include "ChiralViolationContrib.h"
7 #include "ChiralSet.h"
8 #include <ForceField/ForceField.h>
9
10 namespace DistGeom {
ChiralViolationContrib(ForceFields::ForceField * owner,const ChiralSet * cset,double weight)11 ChiralViolationContrib::ChiralViolationContrib(ForceFields::ForceField *owner,
12 const ChiralSet *cset,
13 double weight) {
14 PRECONDITION(owner, "bad owner");
15 PRECONDITION(cset, "bad chiral set")
16
17 URANGE_CHECK(cset->d_idx1, owner->positions().size());
18 URANGE_CHECK(cset->d_idx2, owner->positions().size());
19 URANGE_CHECK(cset->d_idx3, owner->positions().size());
20 URANGE_CHECK(cset->d_idx4, owner->positions().size());
21
22 dp_forceField = owner;
23
24 d_idx1 = cset->d_idx1;
25 d_idx2 = cset->d_idx2;
26 d_idx3 = cset->d_idx3;
27 d_idx4 = cset->d_idx4;
28
29 d_volLower = cset->getLowerVolumeBound();
30 d_volUpper = cset->getUpperVolumeBound();
31
32 d_weight = weight;
33 }
34
getEnergy(double * pos) const35 double ChiralViolationContrib::getEnergy(double *pos) const {
36 PRECONDITION(dp_forceField, "no owner");
37 PRECONDITION(pos, "bad vector");
38
39 unsigned int dim = dp_forceField->dimension();
40 double vol = calcChiralVolume(d_idx1, d_idx2, d_idx3, d_idx4, pos, dim);
41 double res = 0.0;
42 if (vol < d_volLower) {
43 res = d_weight * (vol - d_volLower) * (vol - d_volLower);
44 } else if (vol > d_volUpper) {
45 res = d_weight * (vol - d_volUpper) * (vol - d_volUpper);
46 }
47 // std::cerr<<"Chiral Violation vol: "<<vol<<" E: "<<res<<std::endl;
48 return res;
49 }
50
getGrad(double * pos,double * grad) const51 void ChiralViolationContrib::getGrad(double *pos, double *grad) const {
52 PRECONDITION(dp_forceField, "no owner");
53 PRECONDITION(pos, "bad vector");
54
55 unsigned int dim = dp_forceField->dimension();
56
57 // even if we are minimizing in higher dimension the chiral volume is
58 // calculated using only the first 3 dimensions
59 RDGeom::Point3D v1(pos[d_idx1 * dim] - pos[d_idx4 * dim],
60 pos[d_idx1 * dim + 1] - pos[d_idx4 * dim + 1],
61 pos[d_idx1 * dim + 2] - pos[d_idx4 * dim + 2]);
62
63 RDGeom::Point3D v2(pos[d_idx2 * dim] - pos[d_idx4 * dim],
64 pos[d_idx2 * dim + 1] - pos[d_idx4 * dim + 1],
65 pos[d_idx2 * dim + 2] - pos[d_idx4 * dim + 2]);
66
67 RDGeom::Point3D v3(pos[d_idx3 * dim] - pos[d_idx4 * dim],
68 pos[d_idx3 * dim + 1] - pos[d_idx4 * dim + 1],
69 pos[d_idx3 * dim + 2] - pos[d_idx4 * dim + 2]);
70
71 RDGeom::Point3D v2xv3 = v2.crossProduct(v3);
72
73 double vol = v1.dotProduct(v2xv3);
74 double preFactor;
75
76 if (vol < d_volLower) {
77 preFactor = d_weight * (vol - d_volLower);
78 } else if (vol > d_volUpper) {
79 preFactor = d_weight * (vol - d_volUpper);
80 } else {
81 return;
82 }
83
84 // now comes the hard part - there are a total of 12 variables involved
85 // 4 x 3 - four points and 3 dimensions
86 //
87 grad[dim * d_idx1] += preFactor * ((v2.y) * (v3.z) - (v3.y) * (v2.z));
88 grad[dim * d_idx1 + 1] += preFactor * ((v3.x) * (v2.z) - (v2.x) * (v3.z));
89 grad[dim * d_idx1 + 2] += preFactor * ((v2.x) * (v3.y) - (v3.x) * (v2.y));
90
91 grad[dim * d_idx2] += preFactor * ((v3.y) * (v1.z) - (v3.z) * (v1.y));
92 grad[dim * d_idx2 + 1] += preFactor * ((v3.z) * (v1.x) - (v3.x) * (v1.z));
93 grad[dim * d_idx2 + 2] += preFactor * ((v3.x) * (v1.y) - (v3.y) * (v1.x));
94
95 grad[dim * d_idx3] += preFactor * ((v2.z) * (v1.y) - (v2.y) * (v1.z));
96 grad[dim * d_idx3 + 1] += preFactor * ((v2.x) * (v1.z) - (v2.z) * (v1.x));
97 grad[dim * d_idx3 + 2] += preFactor * ((v2.y) * (v1.x) - (v2.x) * (v1.y));
98
99 grad[dim * d_idx4] +=
100 preFactor *
101 (pos[d_idx1 * dim + 2] * (pos[d_idx2 * dim + 1] - pos[d_idx3 * dim + 1]) +
102 pos[d_idx2 * dim + 2] * (pos[d_idx3 * dim + 1] - pos[d_idx1 * dim + 1]) +
103 pos[d_idx3 * dim + 2] * (pos[d_idx1 * dim + 1] - pos[d_idx2 * dim + 1]));
104
105 grad[dim * d_idx4 + 1] +=
106 preFactor *
107 (pos[d_idx1 * dim] * (pos[d_idx2 * dim + 2] - pos[d_idx3 * dim + 2]) +
108 pos[d_idx2 * dim] * (pos[d_idx3 * dim + 2] - pos[d_idx1 * dim + 2]) +
109 pos[d_idx3 * dim] * (pos[d_idx1 * dim + 2] - pos[d_idx2 * dim + 2]));
110
111 grad[dim * d_idx4 + 2] +=
112 preFactor *
113 (pos[d_idx1 * dim + 1] * (pos[d_idx2 * dim] - pos[d_idx3 * dim]) +
114 pos[d_idx2 * dim + 1] * (pos[d_idx3 * dim] - pos[d_idx1 * dim]) +
115 pos[d_idx3 * dim + 1] * (pos[d_idx1 * dim] - pos[d_idx2 * dim]));
116 // std::cerr<<"Chiral Violation grad: "<<preFactor<<std::endl;
117 }
118 }
119