-
Notifications
You must be signed in to change notification settings - Fork 0
/
bonddihedtopology.C
217 lines (161 loc) · 9.49 KB
/
bonddihedtopology.C
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
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
#include <stdio.h>
#include <iostream>
#include <math.h>
#include <stdlib.h>
#include <string.h>
#include <fstream>
#include <sstream>
#include <string>
#include "FreeEnergyHS.h"
#include "Tools.h"
using namespace std;
void hsc::bond_old_topology(){
if ( open_chain ){
if ( clusterEnd_1 <= clusterStart_1 ){
oldEBond_1 = bond_bprotateEnergy(&part[clusterStart_1], &part[(N/2) - 1]);
oldEBond_1 += bond_bprotateEnergy(&part[0], &part[clusterEnd_1 + 1]);
}else{
oldEBond_1 = bond_bprotateEnergy(&part[clusterStart_1], &part[clusterEnd_1]);
}
if ( clusterEnd_2 <= clusterStart_2 ){
oldEBond_2 = bond_bprotateEnergy(&part[clusterStart_2], &part[N-1]);
oldEBond_2 += bond_bprotateEnergy(&part[N/2], &part[clusterEnd_2 + 1]);
}else{
oldEBond_2 = bond_bprotateEnergy(&part[clusterEnd_2], &part[clusterStart_2 + 1]);
}
}else{
if ( clusterEnd_1 < clusterStart_1 && clusterEnd_2 < clusterStart_2 && clusterEnd_1 == 0 ){
oldEBond_1 = bond_bprotateEnergy(&part[clusterStart_1 - 2], &part[(N/2) - 2]);
oldEBond_2 = bond_bprotateEnergy(&part[N - 2], &part[clusterEnd_2 + 2]);
oldEBond_2 += bond_bprotateEnergy(&part[clusterStart_2], &part[N - 2]);
} else if ( clusterEnd_1 < clusterStart_1 && clusterStart_2 < clusterEnd_2){
oldEBond_1 = bond_bprotateEnergy(&part[clusterStart_1 - 2], &part[0]);
oldEBond_1 += bond_bprotateEnergy(&part[0], &part[clusterEnd_1]);
oldEBond_2 = bond_bprotateEnergy(&part[clusterEnd_2], &part[N - 2]);
oldEBond_2 += bond_bprotateEnergy(&part[N - 2], &part[clusterStart_2 + 2]);
}else if (clusterStart_1 == 0 ){
oldEBond_1 = bond_bprotateEnergy(&part[clusterStart_1], &part[clusterEnd_1]);
oldEBond_2 = bond_bprotateEnergy(&part[clusterStart_2], &part[clusterEnd_2]);
}else if ( clusterStart_2 < clusterEnd_2 ) {
oldEBond_1 = bond_bprotateEnergy(&part[clusterStart_1 - 2], &part[clusterEnd_1]);
oldEBond_2 = bond_bprotateEnergy(&part[clusterStart_2], &part[N - 2]);
oldEBond_2 += bond_bprotateEnergy(&part[N/2], &part[clusterEnd_2]);
} else {
oldEBond_1 = bond_bprotateEnergy(&part[clusterStart_1 - 2], &part[clusterEnd_1]);
oldEBond_2 = bond_bprotateEnergy(&part[clusterEnd_2], &part[clusterStart_2 + 2]);
}
}
}
void hsc::bond_new_topology(){
if ( open_chain ){
if ( clusterEnd_1 <= clusterStart_1 ){
newEBond_1 = bond_bprotateEnergy(&part[clusterStart_1], &part[(N/2) - 1]);
newEBond_1 += bond_bprotateEnergy(&part[0], &part[clusterEnd_1 + 1]);
}else{
newEBond_1 = bond_bprotateEnergy(&part[clusterStart_1], &part[clusterEnd_1]);
}
if ( clusterEnd_2 <= clusterStart_2 ){
newEBond_2 = bond_bprotateEnergy(&part[clusterStart_2], &part[N-1]);
newEBond_2 += bond_bprotateEnergy(&part[N/2], &part[clusterEnd_2 + 1]);
}else{
newEBond_2 = bond_bprotateEnergy(&part[clusterEnd_2], &part[clusterStart_2 + 1]);
}
}else{
if ( clusterEnd_1 < clusterStart_1 && clusterEnd_2 < clusterStart_2 && clusterEnd_1 == 0 ){
newEBond_1 = bond_bprotateEnergy(&part[clusterStart_1 - 2], &part[(N/2) - 2]);
newEBond_2 = bond_bprotateEnergy(&part[N - 2], &part[clusterEnd_2 + 2]);
newEBond_2 += bond_bprotateEnergy(&part[clusterStart_2], &part[N - 2]);
} else if ( clusterEnd_1 < clusterStart_1 && clusterStart_2 < clusterEnd_2){
newEBond_1 = bond_bprotateEnergy(&part[clusterStart_1 - 2], &part[0]);
newEBond_1 += bond_bprotateEnergy(&part[0], &part[clusterEnd_1]);
newEBond_2 = bond_bprotateEnergy(&part[clusterEnd_2], &part[N - 2]);
newEBond_2 += bond_bprotateEnergy(&part[N - 2], &part[clusterStart_2 + 2]);
}else if (clusterStart_1 == 0 ){
newEBond_1 = bond_bprotateEnergy(&part[clusterStart_1], &part[clusterEnd_1]);
newEBond_2 = bond_bprotateEnergy(&part[clusterStart_2], &part[clusterEnd_2]);
}else if ( clusterStart_2 < clusterEnd_2 ) {
newEBond_1 = bond_bprotateEnergy(&part[clusterStart_1 - 2], &part[clusterEnd_1]);
newEBond_2 = bond_bprotateEnergy(&part[clusterStart_2], &part[N - 2]);
newEBond_2 += bond_bprotateEnergy(&part[N/2], &part[clusterEnd_2]);
} else {
newEBond_1 = bond_bprotateEnergy(&part[clusterStart_1 - 2], &part[clusterEnd_1]);
newEBond_2 = bond_bprotateEnergy(&part[clusterEnd_2], &part[clusterStart_2 + 2]);
}
}
}
void hsc::dihed_old_topology(){
if ( open_chain ) {
if( clusterEnd_1 <= clusterStart_1 ){ //account for a break in the chain.
oldeDihed_1 = dihed_forward(&part[clusterStart_1], &part[(N/2) - 1]);
oldeDihed_1 += dihed_forward(&part[0], &part[clusterEnd_1 + 1]);
}else{
oldeDihed_1 = dihed_forward(&part[clusterStart_1], &part[clusterEnd_1]);
}
if( clusterEnd_2 <= clusterStart_2 ){ //account for a break in the chain.
oldeDihed_2 = dihed_forward(&part[clusterStart_2], &part[N-1]);
oldeDihed_2 += dihed_forward(&part[N/2], &part[clusterEnd_2 + 1]);
}else{
oldeDihed_2 = dihed_forward(&part[clusterStart_2], &part[clusterEnd_2 + 1]);
}
}else{
if ( clusterEnd_1 < clusterStart_1 && clusterEnd_2 < clusterStart_2 && clusterEnd_1 == 0 ){
oldeDihed_1 = dihed_forward(&part[clusterStart_1], &part[N/2 - 1]);
oldeDihed_2 = dihed_forward(&part[clusterStart_2 + 2], &part[clusterEnd_2 + 1]);
}
else if ( clusterEnd_1 < clusterStart_1 && clusterStart_2 < clusterEnd_2 ){
oldeDihed_1 = dihed_forward(&part[clusterStart_1], &part[(N/2) - 1]);
oldeDihed_1 += dihed_forward(&part[3], &part[clusterEnd_1 + 1]);
oldeDihed_2 = dihed_forward(&part[clusterEnd_2], &part[(N - 1)]);
oldeDihed_2 += dihed_forward(&part[(N/2) + 3], &part[clusterStart_2 + 1]);
}else if ( clusterEnd_2 == N/2){
oldeDihed_1 = dihed_forward(&part[clusterStart_1], &part[clusterEnd_1 + 1]);
oldeDihed_2 = dihed_forward(&part[clusterEnd_2], &part[(N - 1)]);
oldeDihed_2 += dihed_forward(&part[(N/2 + 3)], &part[clusterStart_2 + 1]);
}else if (clusterStart_1 == 0 ){
oldeDihed_1 = dihed_forward(&part[clusterStart_1], &part[N/2 - 1]);
oldeDihed_1 += dihed_forward(&part[0], &part[clusterEnd_1 + 1]);
oldeDihed_2 = dihed_forward(&part[clusterStart_2], &part[clusterEnd_2 + 1]);
}else{
oldeDihed_1 = dihed_forward(&part[clusterStart_1], &part[clusterEnd_1 + 1]);
oldeDihed_2 = dihed_forward(&part[clusterEnd_2], &part[clusterStart_2 + 1]);
}
}
}
void hsc::dihed_new_topolgy(){
if ( open_chain ) {
if( clusterEnd_1 <= clusterStart_1 ){ //account for a break in the chain.
neweDihed_1 = dihed_forward(&part[clusterStart_1], &part[(N/2) - 1]);
neweDihed_1 += dihed_forward(&part[0], &part[clusterEnd_1 + 1]);
}else{
neweDihed_1 = dihed_forward(&part[clusterStart_1], &part[clusterEnd_1]);
}
if( clusterEnd_2 <= clusterStart_2 ){ //account for a break in the chain.
neweDihed_2 = dihed_forward(&part[clusterStart_2], &part[N-1]);
neweDihed_2 += dihed_forward(&part[N/2], &part[clusterEnd_2 + 1]);
}else{
neweDihed_2 = dihed_forward(&part[clusterStart_2], &part[clusterEnd_2 + 1]);
}
}else{
if ( clusterEnd_1 < clusterStart_1 && clusterEnd_2 < clusterStart_2 && clusterEnd_1 == 0 ){
neweDihed_1 = dihed_forward(&part[clusterStart_1], &part[N/2 - 1]);
neweDihed_2 = dihed_forward(&part[clusterStart_2 + 2], &part[clusterEnd_2 + 1]);
}
else if ( clusterEnd_1 < clusterStart_1 && clusterStart_2 < clusterEnd_2 ){
neweDihed_1 = dihed_forward(&part[clusterStart_1], &part[(N/2) - 1]);
neweDihed_1 += dihed_forward(&part[3], &part[clusterEnd_1 + 1]);
neweDihed_2 = dihed_forward(&part[clusterEnd_2], &part[(N - 1)]);
neweDihed_2 += dihed_forward(&part[(N/2) + 3], &part[clusterStart_2 + 1]);
}else if ( clusterEnd_2 == N/2){
neweDihed_1 = dihed_forward(&part[clusterStart_1], &part[clusterEnd_1 + 1]);
neweDihed_2 = dihed_forward(&part[clusterEnd_2], &part[N - 1]);
neweDihed_2 += dihed_forward(&part[N/2 + 3], &part[clusterStart_2 + 1]);
}else if (clusterStart_1 == 0 ){
neweDihed_1 = dihed_forward(&part[clusterStart_1], &part[N/2 - 1]);
neweDihed_2 = dihed_forward(&part[clusterEnd_2], &part[(N - 1)]);
neweDihed_2 += dihed_forward(&part[(N/2)], &part[clusterStart_2 + 1]);
}else{
neweDihed_1 = dihed_forward(&part[clusterStart_1], &part[clusterEnd_1 + 1]);
neweDihed_2 = dihed_forward(&part[clusterEnd_2], &part[clusterStart_2 + 1]);
}
}
}