-
Notifications
You must be signed in to change notification settings - Fork 0
/
Copy pathabfe_restraints_nb.txt
101 lines (77 loc) · 3.06 KB
/
abfe_restraints_nb.txt
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
(* I. === Some definitions === *)
kt = 300. *8.314/4184 (* reasonable value for R*T at 300K in kcal/
mol; rerun everything if this input is changed *)
beta = 1/kt
vol0 = 1660. (* standard volume for 1mol/L *)
(* II. === Evaluation of the integrals *)
(* == dihedral terms == *)
(* RRHO *)
c1all = Integrate[Exp[-b k0 s^2], {s, -Infinity, Infinity},
Assumptions -> b > 0 && k0 > 0 && x0 > 0]
(* Schroedinger's expression, after variable substitution*)
c1afullvarsubst =
Integrate[Exp[-b k0 s^2], {s, -Pi, Pi},
Assumptions -> b > 0 && k0 > 0]
(* == distance restraint == *)
(* Schroedinger's expression *)
c2full =
Integrate[x^2 Exp[-b k0 (x - x0)^2], {x, 0, Infinity},
Assumptions -> b > 0 && k0 > 0 ]
(* RRHO *)
c2rr = x0^2 Integrate[ Exp[-b k0 s^2], {s, -Infinity, Infinity},
Assumptions -> b > 0 && k0 > 0]
(* == angle restraint == *)
(* exact -- brute force *)
c3full =
Integrate[Sin[x] Exp[-b k0 (x - x0)^2], {x, 0, Pi},
Assumptions -> b > 0 && k0 > 0 && x0 > 0]
(* Schroedinger's expression, MMA needs variable subst. *)
c3all = Integrate[
Sin[s + x0] Exp[-b k0 s^2], {s, -Infinity, Infinity},
Assumptions -> b > 0 && k0 > 0 && x0 > 0]
(* RRHO *)
c3rr = Sin[x0] Integrate[Exp[-b k0 s^2], {s, -Infinity, Infinity},
Assumptions -> b > 0 && k0 > 0]
(* III. === Numerical examples === *)
(* == L1A-5 == *)
(* RRHO *)
-kt Log[(8 Pi^2 vol0)/((c2rr /. {b -> beta, x0 -> 5.1, k0 -> 5.})*
(c3rr /. {b -> beta, k0 -> 5, x0 -> 67.5 Degree})*
(c3rr /. {b -> beta, k0 -> 5, x0 -> 84.5 Degree})*
(c1all /. {b -> beta, k0 -> 2.5})*
(c1all /. {b -> beta, k0 -> 2.5})*
(c1all /. {b -> beta, k0 -> 2.5})
)]
(* Schroedinger *)
-kt Log[(8 Pi^2 vol0)/((c2full /. {b -> beta, x0 -> 5.1, k0 -> 5.})*
(c3all /. {b -> beta, k0 -> 5., x0 -> 67.5 Degree})*
(c3all /. {b -> beta, k0 -> 5., x0 -> 84.5 Degree})*
(c1afullvarsubst /. {b -> beta, k0 -> 2.5})*
(c1afullvarsubst /. {b -> beta, k0 -> 2.5})*
(c1afullvarsubst /. {b -> beta, k0 -> 2.5})
)]
(* == The EXTR example == *)
(* RRHO *)
-kt Log[(8 Pi^2 vol0)/((c2rr /. {b -> beta, x0 -> 3., k0 -> 1.})*
(c3rr /. {b -> beta, k0 -> 1., x0 -> 22.5 Degree})*
(c3rr /. {b -> beta, k0 -> 1., x0 -> (180 - 22.5) Degree})*
(c1all /. {b -> beta, k0 -> 1.})*
(c1all /. {b -> beta, k0 -> 1.})*
(c1all /. {b -> beta, k0 -> 1.})
)]
(* Schroedinger *)
-kt Log[(8 Pi^2 vol0)/((c2full /. {b -> beta, x0 -> 3., k0 -> 1.})*
(c3all /. {b -> beta, k0 -> 1., x0 -> 22.5 Degree})*
(c3all /. {b -> beta, k0 -> 1., x0 -> (180 - 22.5) Degree})*
(c1afullvarsubst /. {b -> beta, k0 -> 1.})*
(c1afullvarsubst /. {b -> beta, k0 -> 1.})*
(c1afullvarsubst /. {b -> beta, k0 -> 1.})
)]
(* and fully exact *)
-kt Log[(8 Pi^2 vol0)/((c2full /. {b -> beta, x0 -> 3., k0 -> 1.})*
(Re[c3full /. {b -> beta, k0 -> 1., x0 -> 22.5 Degree}])*
(Re[c3full /. {b -> beta, k0 -> 1., x0 -> (180 - 22.5) Degree}])*
(c1afullvarsubst /. {b -> beta, k0 -> 1.})*
(c1afullvarsubst /. {b -> beta, k0 -> 1.})*
(c1afullvarsubst /. {b -> beta, k0 -> 1.})
)]