Skip to content

Commit 37c2a55

Browse files
drmeisterdrmeister
committed
Added two more examples
These are only partially implemented.
1 parent 5e7571b commit 37c2a55

27 files changed

+6985
-0
lines changed

‎tlm3_dihedral/_DihedralEnergy_debugEvalDeclares.cc

Lines changed: 22 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,22 @@
1+
// DEFINE_EVAL
2+
struct {
3+
double x1;
4+
double y1;
5+
double z1;
6+
double x2;
7+
double y2;
8+
double z2;
9+
double x3;
10+
double y3;
11+
double z3;
12+
double x4;
13+
double y4;
14+
double z4;
15+
double V;
16+
double DN;
17+
double IN;
18+
double cosPhase;
19+
double sinPhase;
20+
double Energy;
21+
double DihedralDeviation;
22+
} eval;

‎tlm3_dihedral/_DihedralEnergy_debugEvalSerialize.cc

Lines changed: 21 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,21 @@
1+
// XML_EVAL
2+
node->attributeIfNotDefault("evalx1", this->eval.x1, 0.0 );
3+
node->attributeIfNotDefault("evaly1", this->eval.y1, 0.0 );
4+
node->attributeIfNotDefault("evalz1", this->eval.z1, 0.0 );
5+
node->attributeIfNotDefault("evalx2", this->eval.x2, 0.0 );
6+
node->attributeIfNotDefault("evaly2", this->eval.y2, 0.0 );
7+
node->attributeIfNotDefault("evalz2", this->eval.z2, 0.0 );
8+
node->attributeIfNotDefault("evalx3", this->eval.x3, 0.0 );
9+
node->attributeIfNotDefault("evaly3", this->eval.y3, 0.0 );
10+
node->attributeIfNotDefault("evalz3", this->eval.z3, 0.0 );
11+
node->attributeIfNotDefault("evalx4", this->eval.x4, 0.0 );
12+
node->attributeIfNotDefault("evaly4", this->eval.y4, 0.0 );
13+
node->attributeIfNotDefault("evalz4", this->eval.z4, 0.0 );
14+
node->attributeIfNotDefault("evalV", this->eval.V, 0.0 );
15+
node->attributeIfNotDefault("evalDN", this->eval.DN, 0.0 );
16+
node->attributeIfNotDefault("evalIN", this->eval.IN, 0.0 );
17+
node->attributeIfNotDefault("evalcosPhase", this->eval.cosPhase, 0.0 );
18+
node->attributeIfNotDefault("evalsinPhase", this->eval.sinPhase, 0.0 );
19+
node->attributeIfNotDefault("evalEnergy", this->eval.Energy, 0.0 );
20+
node->attributeIfNotDefault("evalDihedralDeviation", this->eval.DihedralDeviation, 0.0 );
21+

‎tlm3_dihedral/_DihedralEnergy_debugEvalSet.cc

Lines changed: 21 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,21 @@
1+
// SET_EVAL
2+
EVAL_SET(x1, x1);
3+
EVAL_SET(y1, y1);
4+
EVAL_SET(z1, z1);
5+
EVAL_SET(x2, x2);
6+
EVAL_SET(y2, y2);
7+
EVAL_SET(z2, z2);
8+
EVAL_SET(x3, x3);
9+
EVAL_SET(y3, y3);
10+
EVAL_SET(z3, z3);
11+
EVAL_SET(x4, x4);
12+
EVAL_SET(y4, y4);
13+
EVAL_SET(z4, z4);
14+
EVAL_SET(V, V);
15+
EVAL_SET(DN, DN);
16+
EVAL_SET(IN, IN);
17+
EVAL_SET(cosPhase, cosPhase);
18+
EVAL_SET(sinPhase, sinPhase);
19+
EVAL_SET(Energy, Energy);
20+
EVAL_SET(DihedralDeviation, DihedralDeviation);
21+

‎tlm3_dihedral/_DihedralEnergy_debugFiniteDifference.cc

Lines changed: 661 additions & 0 deletions
Large diffs are not rendered by default.

‎tlm3_dihedral/_DihedralEnergy_math.m

Lines changed: 141 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,141 @@
1+
{Name -> "DihedralEnergy", AdditionalCDeclares -> "",
2+
Input -> {x1, y1, z1, x2, y2, z2, x3, y3, z3, x4, y4, z4, V, DN, IN,
3+
cosPhase, sinPhase}, Output -> {Energy, DihedralDeviation},
4+
EnergyFunction -> "NotUsed", DerivativeVariables ->
5+
{x1, y1, z1, x2, y2, z2, x3, y3, z3, x4, y4, z4},
6+
Rules -> {CCode["DIHEDRAL_SET_PARAMETER(sinPhase);"],
7+
CCode["DIHEDRAL_SET_PARAMETER(cosPhase);"],
8+
CCode["DIHEDRAL_SET_PARAMETER(V);"], CCode["DIHEDRAL_SET_PARAMETER(DN);"],
9+
CCode["DIHEDRAL_SET_PARAMETER(IN);"],
10+
CCode["DIHEDRAL_SET_PARAMETER(I1);"],
11+
CCode["DIHEDRAL_SET_PARAMETER(I2);"],
12+
CCode["DIHEDRAL_SET_PARAMETER(I3);"],
13+
CCode["DIHEDRAL_SET_PARAMETER(I4);"],
14+
CCode["DIHEDRAL_SET_POSITION(x1,I1,0);"],
15+
CCode["DIHEDRAL_SET_POSITION(y1,I1,1);"],
16+
CCode["DIHEDRAL_SET_POSITION(z1,I1,2);"],
17+
CCode["DIHEDRAL_SET_POSITION(x2,I2,0);"],
18+
CCode["DIHEDRAL_SET_POSITION(y2,I2,1);"],
19+
CCode["DIHEDRAL_SET_POSITION(z2,I2,2);"],
20+
CCode["DIHEDRAL_SET_POSITION(x3,I3,0);"],
21+
CCode["DIHEDRAL_SET_POSITION(y3,I3,1);"],
22+
CCode["DIHEDRAL_SET_POSITION(z3,I3,2);"],
23+
CCode["DIHEDRAL_SET_POSITION(x4,I4,0);"],
24+
CCode["DIHEDRAL_SET_POSITION(y4,I4,1);"],
25+
CCode["DIHEDRAL_SET_POSITION(z4,I4,2);"], -(x2*y1) -> tx69, x3*y1 -> tx70,
26+
x1*y2 -> tx71, -(x3*y2) -> tx72, -(x1*y3) -> tx73, x2*y3 -> tx74,
27+
x2*z1 -> tx75, -(x3*z1) -> tx76, -(y2*z1) -> tx77, y3*z1 -> tx78,
28+
-(x1*z2) -> tx79, x3*z2 -> tx80, y1*z2 -> tx81, -(y3*z2) -> tx82,
29+
x1*z3 -> tx83, -(x2*z3) -> tx84, -(y1*z3) -> tx85, y2*z3 -> tx86,
30+
tx69 + tx70 + tx71 + tx72 + tx73 + tx74 -> tx87,
31+
tx75 + tx76 + tx79 + tx80 + tx83 + tx84 -> tx88,
32+
tx77 + tx78 + tx81 + tx82 + tx85 + tx86 -> tx89, power2[tx87] -> tx90,
33+
power2[tx88] -> tx91, power2[tx89] -> tx92, tx90 + tx91 + tx92 -> tx93,
34+
mysqrt[tx93] -> LenA, x4*y2 -> tx94, -(x4*y3) -> tx95, -(x2*y4) -> tx96,
35+
x3*y4 -> tx97, -(x4*z2) -> tx98, y4*z2 -> tx99, x4*z3 -> tx100,
36+
-(y4*z3) -> tx101, x2*z4 -> tx102, -(x3*z4) -> tx103, -(y2*z4) -> tx104,
37+
y3*z4 -> tx105, tx72 + tx74 + tx94 + tx95 + tx96 + tx97 -> tx106,
38+
tx100 + tx102 + tx103 + tx80 + tx84 + tx98 -> tx107,
39+
tx101 + tx104 + tx105 + tx82 + tx86 + tx99 -> tx108,
40+
power2[tx106] -> tx109, power2[tx107] -> tx110, power2[tx108] -> tx111,
41+
tx109 + tx110 + tx111 -> tx112, mysqrt[tx112] -> LenB,
42+
reciprocal[LenA] -> ReciprocalLenA, reciprocal[LenB] -> ReciprocalLenB,
43+
CCode["if (fabs(LenA)<TENM3) ReciprocalLenA = 0.0;"],
44+
CCode["if (fabs(LenB)<TENM3) ReciprocalLenB = 0.0;"],
45+
ReciprocalLenA*ReciprocalLenB -> RecLenARecLenB,
46+
CCode["EraseLinearDihedral = 1.0;"],
47+
CCode["if (RecLenARecLenB==0.0) EraseLinearDihedral = 0.0;"],
48+
tx106*tx87 -> tx113, tx107*tx88 -> tx114, tx108*tx89 -> tx115,
49+
tx113 + tx114 + tx115 -> tx116, RecLenARecLenB*tx116 -> CosPhi,
50+
-x3 -> tx117, -y3 -> tx118, -z3 -> tx119, tx117 + x2 -> tx120,
51+
tx118 + y2 -> tx121, tx119 + z2 -> tx122, power2[tx120] -> tx123,
52+
power2[tx121] -> tx124, power2[tx122] -> tx125, tx117 + x4 -> tx126,
53+
tx118 + y4 -> tx127, tx119 + z4 -> tx128, tx123 + tx124 + tx125 -> tx129,
54+
tx128*tx87 -> tx130, tx127*tx88 -> tx131, tx126*tx89 -> tx132,
55+
mysqrt[tx129] -> tx133, tx130 + tx131 + tx132 -> tx134,
56+
RecLenARecLenB*tx133*tx134 -> SinPhi,
57+
CCode["CosPhi=MAX(-1.0,MIN(1.0,CosPhi));"],
58+
MathCode["CosNPhi = mathCosNPhi[IN,SinPhi,CosPhi];"],
59+
MathCode["SinNPhi = mathSinNPhi[IN,SinPhi,CosPhi];"],
60+
CCode["sinNPhiCosNPhi(IN, &SinNPhi, &CosNPhi, SinPhi, CosPhi);"],
61+
CosNPhi*cosPhase -> tx135, SinNPhi*sinPhase -> tx136,
62+
1. + tx135 + tx136 -> DihedralDeviation,
63+
DihedralDeviation*EraseLinearDihedral*V -> Energy,
64+
CCode["DIHEDRAL_ENERGY_ACCUMULATE(Energy);"]},
65+
HessianStructure -> {{14, 26, 27, 28, 29, 30, 31, 32, 33, 34, 35, 36},
66+
{26, 15, 37, 38, 39, 40, 41, 42, 43, 44, 45, 46},
67+
{27, 37, 16, 47, 48, 49, 50, 51, 52, 53, 54, 55},
68+
{28, 38, 47, 17, 56, 57, 58, 59, 60, 61, 62, 63},
69+
{29, 39, 48, 56, 18, 64, 65, 66, 67, 68, 69, 70},
70+
{30, 40, 49, 57, 64, 19, 71, 72, 73, 74, 75, 76},
71+
{31, 41, 50, 58, 65, 71, 20, 77, 78, 79, 80, 81},
72+
{32, 42, 51, 59, 66, 72, 77, 21, 82, 83, 84, 85},
73+
{33, 43, 52, 60, 67, 73, 78, 82, 22, 86, 87, 88},
74+
{34, 44, 53, 61, 68, 74, 79, 83, 86, 23, 89, 90},
75+
{35, 45, 54, 62, 69, 75, 80, 84, 87, 89, 24, 91},
76+
{36, 46, 55, 63, 70, 76, 81, 85, 88, 90, 91, 25}}, CCode -> "NOT USED",
77+
MathCode -> "Block[{CosPhi,DihedralDeviation,Energy,LenA,LenB,ReciprocalLenA\
78+
,ReciprocalLenB,RecLenARecLenB,SinPhi,tx100,tx101,tx102,tx103,tx104,tx105,tx1\
79+
06,tx107,tx108,tx109,tx110,tx111,tx112,tx113,tx114,tx115,tx116,tx117,tx118,tx\
80+
119,tx120,tx121,tx122,tx123,tx124,tx125,tx126,tx127,tx128,tx129,tx130,tx131,t\
81+
x132,tx133,tx134,tx135,tx136,tx69,tx70,tx71,tx72,tx73,tx74,tx75,tx76,tx77,tx7\
82+
8,tx79,tx80,tx81,tx82,tx83,tx84,tx85,tx86,tx87,tx88,tx89,tx90,tx91,tx92,tx93,\
83+
tx94,tx95,tx96,tx97,tx98,tx99,xxxDummy},\n\t(*DIHEDRAL_SET_PARAMETER(sinPhase\
84+
);*)\n\t(*DIHEDRAL_SET_PARAMETER(cosPhase);*)\n\t(*DIHEDRAL_SET_PARAMETER(V);\
85+
*)\n\t(*DIHEDRAL_SET_PARAMETER(DN);*)\n\t(*DIHEDRAL_SET_PARAMETER(IN);*)\n\t(\
86+
*DIHEDRAL_SET_PARAMETER(I1);*)\n\t(*DIHEDRAL_SET_PARAMETER(I2);*)\n\t(*DIHEDR\
87+
AL_SET_PARAMETER(I3);*)\n\t(*DIHEDRAL_SET_PARAMETER(I4);*)\n\t(*DIHEDRAL_SET_\
88+
POSITION(x1,I1,0);*)\n\t(*DIHEDRAL_SET_POSITION(y1,I1,1);*)\n\t(*DIHEDRAL_SET\
89+
_POSITION(z1,I1,2);*)\n\t(*DIHEDRAL_SET_POSITION(x2,I2,0);*)\n\t(*DIHEDRAL_SE\
90+
T_POSITION(y2,I2,1);*)\n\t(*DIHEDRAL_SET_POSITION(z2,I2,2);*)\n\t(*DIHEDRAL_S\
91+
ET_POSITION(x3,I3,0);*)\n\t(*DIHEDRAL_SET_POSITION(y3,I3,1);*)\n\t(*DIHEDRAL_\
92+
SET_POSITION(z3,I3,2);*)\n\t(*DIHEDRAL_SET_POSITION(x4,I4,0);*)\n\t(*DIHEDRAL\
93+
_SET_POSITION(y4,I4,1);*)\n\t(*DIHEDRAL_SET_POSITION(z4,I4,2);*)\n\ttx69=-(x2 \
94+
y1); (*Rule 22*)\n\ttx70=x3 y1; (*Rule 23*)\n\ttx71=x1 y2; (*Rule \
95+
24*)\n\ttx72=-(x3 y2); (*Rule 25*)\n\ttx73=-(x1 y3); (*Rule 26*)\n\ttx74=x2 \
96+
y3; (*Rule 27*)\n\ttx75=x2 z1; (*Rule 28*)\n\ttx76=-(x3 z1); (*Rule \
97+
29*)\n\ttx77=-(y2 z1); (*Rule 30*)\n\ttx78=y3 z1; (*Rule 31*)\n\ttx79=-(x1 \
98+
z2); (*Rule 32*)\n\ttx80=x3 z2; (*Rule 33*)\n\ttx81=y1 z2; (*Rule \
99+
34*)\n\ttx82=-(y3 z2); (*Rule 35*)\n\ttx83=x1 z3; (*Rule 36*)\n\ttx84=-(x2 \
100+
z3); (*Rule 37*)\n\ttx85=-(y1 z3); (*Rule 38*)\n\ttx86=y2 z3; (*Rule \
101+
39*)\n\ttx87=tx69 + tx70 + tx71 + tx72 + tx73 + tx74; (*Rule \
102+
40*)\n\ttx88=tx75 + tx76 + tx79 + tx80 + tx83 + tx84; (*Rule \
103+
41*)\n\ttx89=tx77 + tx78 + tx81 + tx82 + tx85 + tx86; (*Rule 42*)\n\ttx90= \
104+
mathPower2[tx87]; (*Rule 43*)\n\ttx91= mathPower2[tx88]; (*Rule 44*)\n\ttx92= \
105+
mathPower2[tx89]; (*Rule 45*)\n\ttx93=tx90 + tx91 + tx92; (*Rule \
106+
46*)\n\tLenA= mathSqrt[tx93]; (*Rule 47*)\n\ttx94=x4 y2; (*Rule \
107+
48*)\n\ttx95=-(x4 y3); (*Rule 49*)\n\ttx96=-(x2 y4); (*Rule 50*)\n\ttx97=x3 \
108+
y4; (*Rule 51*)\n\ttx98=-(x4 z2); (*Rule 52*)\n\ttx99=y4 z2; (*Rule \
109+
53*)\n\ttx100=x4 z3; (*Rule 54*)\n\ttx101=-(y4 z3); (*Rule 55*)\n\ttx102=x2 \
110+
z4; (*Rule 56*)\n\ttx103=-(x3 z4); (*Rule 57*)\n\ttx104=-(y2 z4); (*Rule \
111+
58*)\n\ttx105=y3 z4; (*Rule 59*)\n\ttx106=tx72 + tx74 + tx94 + tx95 + tx96 + \
112+
tx97; (*Rule 60*)\n\ttx107=tx100 + tx102 + tx103 + tx80 + tx84 + tx98; (*Rule \
113+
61*)\n\ttx108=tx101 + tx104 + tx105 + tx82 + tx86 + tx99; (*Rule \
114+
62*)\n\ttx109= mathPower2[tx106]; (*Rule 63*)\n\ttx110= mathPower2[tx107]; \
115+
(*Rule 64*)\n\ttx111= mathPower2[tx108]; (*Rule 65*)\n\ttx112=tx109 + tx110 + \
116+
tx111; (*Rule 66*)\n\tLenB= mathSqrt[tx112]; (*Rule 67*)\n\tReciprocalLenA= \
117+
mathReciprocal[LenA]; (*Rule 68*)\n\tReciprocalLenB= mathReciprocal[LenB]; \
118+
(*Rule 69*)\n\t(*if (fabs(LenA)<TENM3) ReciprocalLenA = 0.0;*)\n\t(*if \
119+
(fabs(LenB)<TENM3) ReciprocalLenB = 0.0;*)\n\tRecLenARecLenB=ReciprocalLenA \
120+
ReciprocalLenB; (*Rule 72*)\n\t(*EraseLinearDihedral = 1.0;*)\n\t(*if \
121+
(RecLenARecLenB==0.0) EraseLinearDihedral = 0.0;*)\n\ttx113=tx106 tx87; \
122+
(*Rule 75*)\n\ttx114=tx107 tx88; (*Rule 76*)\n\ttx115=tx108 tx89; (*Rule \
123+
77*)\n\ttx116=tx113 + tx114 + tx115; (*Rule 78*)\n\tCosPhi=RecLenARecLenB \
124+
tx116; (*Rule 79*)\n\ttx117=-x3; (*Rule 80*)\n\ttx118=-y3; (*Rule \
125+
81*)\n\ttx119=-z3; (*Rule 82*)\n\ttx120=tx117 + x2; (*Rule \
126+
83*)\n\ttx121=tx118 + y2; (*Rule 84*)\n\ttx122=tx119 + z2; (*Rule \
127+
85*)\n\ttx123= mathPower2[tx120]; (*Rule 86*)\n\ttx124= mathPower2[tx121]; \
128+
(*Rule 87*)\n\ttx125= mathPower2[tx122]; (*Rule 88*)\n\ttx126=tx117 + x4; \
129+
(*Rule 89*)\n\ttx127=tx118 + y4; (*Rule 90*)\n\ttx128=tx119 + z4; (*Rule \
130+
91*)\n\ttx129=tx123 + tx124 + tx125; (*Rule 92*)\n\ttx130=tx128 tx87; (*Rule \
131+
93*)\n\ttx131=tx127 tx88; (*Rule 94*)\n\ttx132=tx126 tx89; (*Rule \
132+
95*)\n\ttx133= mathSqrt[tx129]; (*Rule 96*)\n\ttx134=tx130 + tx131 + tx132; \
133+
(*Rule 97*)\n\tSinPhi=RecLenARecLenB tx133 tx134; (*Rule \
134+
98*)\n\t(*CosPhi=MAX(-1.0,MIN(1.0,CosPhi));*)\n\tCosNPhi = \
135+
mathCosNPhi[IN,SinPhi,CosPhi];\n\tSinNPhi = \
136+
mathSinNPhi[IN,SinPhi,CosPhi];\n\t(*sinNPhiCosNPhi(IN, &SinNPhi, &CosNPhi, \
137+
SinPhi, CosPhi);*)\n\ttx135=CosNPhi cosPhase; (*Rule 103*)\n\ttx136=SinNPhi \
138+
sinPhase; (*Rule 104*)\n\tDihedralDeviation=1. + tx135 + tx136; (*Rule \
139+
105*)\n\tEnergy=DihedralDeviation EraseLinearDihedral V; (*Rule \
140+
106*)\n\t(*DIHEDRAL_ENERGY_ACCUMULATE(Energy);*)\n{Energy, \
141+
DihedralDeviation}]\n"}

‎tlm3_dihedral/_DihedralEnergy_termCode.cc

Lines changed: 108 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,108 @@
1+
// C-code
2+
DIHEDRAL_SET_PARAMETER(sinPhase);
3+
DIHEDRAL_SET_PARAMETER(cosPhase);
4+
DIHEDRAL_SET_PARAMETER(V);
5+
DIHEDRAL_SET_PARAMETER(DN);
6+
DIHEDRAL_SET_PARAMETER(IN);
7+
DIHEDRAL_SET_PARAMETER(I1);
8+
DIHEDRAL_SET_PARAMETER(I2);
9+
DIHEDRAL_SET_PARAMETER(I3);
10+
DIHEDRAL_SET_PARAMETER(I4);
11+
DIHEDRAL_SET_POSITION(x1,I1,0);
12+
DIHEDRAL_SET_POSITION(y1,I1,1);
13+
DIHEDRAL_SET_POSITION(z1,I1,2);
14+
DIHEDRAL_SET_POSITION(x2,I2,0);
15+
DIHEDRAL_SET_POSITION(y2,I2,1);
16+
DIHEDRAL_SET_POSITION(z2,I2,2);
17+
DIHEDRAL_SET_POSITION(x3,I3,0);
18+
DIHEDRAL_SET_POSITION(y3,I3,1);
19+
DIHEDRAL_SET_POSITION(z3,I3,2);
20+
DIHEDRAL_SET_POSITION(x4,I4,0);
21+
DIHEDRAL_SET_POSITION(y4,I4,1);
22+
DIHEDRAL_SET_POSITION(z4,I4,2);
23+
tx69 = -(x2*y1); /* rule 22 */
24+
tx70 = x3*y1; /* rule 23 */
25+
tx71 = x1*y2; /* rule 24 */
26+
tx72 = -(x3*y2); /* rule 25 */
27+
tx73 = -(x1*y3); /* rule 26 */
28+
tx74 = x2*y3; /* rule 27 */
29+
tx75 = x2*z1; /* rule 28 */
30+
tx76 = -(x3*z1); /* rule 29 */
31+
tx77 = -(y2*z1); /* rule 30 */
32+
tx78 = y3*z1; /* rule 31 */
33+
tx79 = -(x1*z2); /* rule 32 */
34+
tx80 = x3*z2; /* rule 33 */
35+
tx81 = y1*z2; /* rule 34 */
36+
tx82 = -(y3*z2); /* rule 35 */
37+
tx83 = x1*z3; /* rule 36 */
38+
tx84 = -(x2*z3); /* rule 37 */
39+
tx85 = -(y1*z3); /* rule 38 */
40+
tx86 = y2*z3; /* rule 39 */
41+
tx87 = tx69 + tx70 + tx71 + tx72 + tx73 + tx74; /* rule 40 */
42+
tx88 = tx75 + tx76 + tx79 + tx80 + tx83 + tx84; /* rule 41 */
43+
tx89 = tx77 + tx78 + tx81 + tx82 + tx85 + tx86; /* rule 42 */
44+
tx90 = power2(tx87); /* rule 43 */
45+
tx91 = power2(tx88); /* rule 44 */
46+
tx92 = power2(tx89); /* rule 45 */
47+
tx93 = tx90 + tx91 + tx92; /* rule 46 */
48+
LenA = mysqrt(tx93); /* rule 47 */
49+
tx94 = x4*y2; /* rule 48 */
50+
tx95 = -(x4*y3); /* rule 49 */
51+
tx96 = -(x2*y4); /* rule 50 */
52+
tx97 = x3*y4; /* rule 51 */
53+
tx98 = -(x4*z2); /* rule 52 */
54+
tx99 = y4*z2; /* rule 53 */
55+
tx100 = x4*z3; /* rule 54 */
56+
tx101 = -(y4*z3); /* rule 55 */
57+
tx102 = x2*z4; /* rule 56 */
58+
tx103 = -(x3*z4); /* rule 57 */
59+
tx104 = -(y2*z4); /* rule 58 */
60+
tx105 = y3*z4; /* rule 59 */
61+
tx106 = tx72 + tx74 + tx94 + tx95 + tx96 + tx97; /* rule 60 */
62+
tx107 = tx100 + tx102 + tx103 + tx80 + tx84 + tx98; /* rule 61 */
63+
tx108 = tx101 + tx104 + tx105 + tx82 + tx86 + tx99; /* rule 62 */
64+
tx109 = power2(tx106); /* rule 63 */
65+
tx110 = power2(tx107); /* rule 64 */
66+
tx111 = power2(tx108); /* rule 65 */
67+
tx112 = tx109 + tx110 + tx111; /* rule 66 */
68+
LenB = mysqrt(tx112); /* rule 67 */
69+
ReciprocalLenA = reciprocal(LenA); /* rule 68 */
70+
ReciprocalLenB = reciprocal(LenB); /* rule 69 */
71+
if (fabs(LenA)<TENM3) ReciprocalLenA = 0.0;
72+
if (fabs(LenB)<TENM3) ReciprocalLenB = 0.0;
73+
RecLenARecLenB = ReciprocalLenA*ReciprocalLenB; /* rule 72 */
74+
EraseLinearDihedral = 1.0;
75+
if (RecLenARecLenB==0.0) EraseLinearDihedral = 0.0;
76+
tx113 = tx106*tx87; /* rule 75 */
77+
tx114 = tx107*tx88; /* rule 76 */
78+
tx115 = tx108*tx89; /* rule 77 */
79+
tx116 = tx113 + tx114 + tx115; /* rule 78 */
80+
CosPhi = RecLenARecLenB*tx116; /* rule 79 */
81+
tx117 = -x3; /* rule 80 */
82+
tx118 = -y3; /* rule 81 */
83+
tx119 = -z3; /* rule 82 */
84+
tx120 = tx117 + x2; /* rule 83 */
85+
tx121 = tx118 + y2; /* rule 84 */
86+
tx122 = tx119 + z2; /* rule 85 */
87+
tx123 = power2(tx120); /* rule 86 */
88+
tx124 = power2(tx121); /* rule 87 */
89+
tx125 = power2(tx122); /* rule 88 */
90+
tx126 = tx117 + x4; /* rule 89 */
91+
tx127 = tx118 + y4; /* rule 90 */
92+
tx128 = tx119 + z4; /* rule 91 */
93+
tx129 = tx123 + tx124 + tx125; /* rule 92 */
94+
tx130 = tx128*tx87; /* rule 93 */
95+
tx131 = tx127*tx88; /* rule 94 */
96+
tx132 = tx126*tx89; /* rule 95 */
97+
tx133 = mysqrt(tx129); /* rule 96 */
98+
tx134 = tx130 + tx131 + tx132; /* rule 97 */
99+
SinPhi = RecLenARecLenB*tx133*tx134; /* rule 98 */
100+
CosPhi=MAX(-1.0,MIN(1.0,CosPhi));
101+
/*CosNPhi = mathCosNPhi[IN,SinPhi,CosPhi];*/
102+
/*SinNPhi = mathSinNPhi[IN,SinPhi,CosPhi];*/
103+
sinNPhiCosNPhi(IN, &SinNPhi, &CosNPhi, SinPhi, CosPhi);
104+
tx135 = CosNPhi*cosPhase; /* rule 103 */
105+
tx136 = SinNPhi*sinPhase; /* rule 104 */
106+
DihedralDeviation = 1. + tx135 + tx136; /* rule 105 */
107+
Energy = DihedralDeviation*EraseLinearDihedral*V; /* rule 106 */
108+
DIHEDRAL_ENERGY_ACCUMULATE(Energy);

0 commit comments

Comments
 (0)