You can not select more than 25 topics Topics must start with a chinese character,a letter or number, can include dashes ('-') and can be up to 35 characters long.

dihedral.py 8.4 kB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186
  1. # Copyright 2021 Huawei Technologies Co., Ltd
  2. #
  3. # Licensed under the Apache License, Version 2.0 (the "License");
  4. # you may not use this file except in compliance with the License.
  5. # You may obtain a copy of the License at
  6. #
  7. # http://www.apache.org/licenses/LICENSE-2.0
  8. #
  9. # Unless required by applicable law or agreed to in writing, software
  10. # distributed under the License is distributed on an "AS IS" BASIS,
  11. # WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
  12. # See the License for the specific language governing permissions and
  13. # limitations under the License.
  14. # ============================================================================
  15. '''Dihedral'''
  16. import math
  17. class Dihedral:
  18. '''Dihedral'''
  19. def __init__(self, controller):
  20. self.CONSTANT_Pi = 3.1415926535897932
  21. if controller.amber_parm is not None:
  22. file_path = controller.amber_parm
  23. self.read_information_from_amberfile(file_path)
  24. def read_information_from_amberfile(self, file_path):
  25. '''read amber file'''
  26. file = open(file_path, 'r')
  27. context = file.readlines()
  28. file.close()
  29. for idx, val in enumerate(context):
  30. if idx < len(context) - 1:
  31. if "%FLAG POINTERS" in val + context[idx + 1] and "%FORMAT(10I8)" in val + context[idx + 1]:
  32. start_idx = idx + 2
  33. count = 0
  34. value = list(map(int, context[start_idx].strip().split()))
  35. self.dihedral_with_hydrogen = value[6]
  36. self.dihedral_numbers = value[7]
  37. self.dihedral_numbers += self.dihedral_with_hydrogen
  38. information = []
  39. information.extend(value)
  40. while count < 15:
  41. start_idx += 1
  42. value = list(map(int, context[start_idx].strip().split()))
  43. information.extend(value)
  44. count += len(value)
  45. self.dihedral_type_numbers = information[17]
  46. print("dihedral type numbers ", self.dihedral_type_numbers)
  47. break
  48. self.phase_type = [0] * self.dihedral_type_numbers
  49. self.pk_type = [0] * self.dihedral_type_numbers
  50. self.pn_type = [0] * self.dihedral_type_numbers
  51. for idx, val in enumerate(context):
  52. if "%FLAG DIHEDRAL_FORCE_CONSTANT" in val:
  53. count = 0
  54. start_idx = idx
  55. information = []
  56. while count < self.dihedral_type_numbers:
  57. start_idx += 1
  58. if "%FORMAT" in context[start_idx]:
  59. continue
  60. else:
  61. value = list(map(float, context[start_idx].strip().split()))
  62. information.extend(value)
  63. count += len(value)
  64. self.pk_type = information[:self.dihedral_type_numbers]
  65. break
  66. for idx, val in enumerate(context):
  67. if "%FLAG DIHEDRAL_PHASE" in val:
  68. count = 0
  69. start_idx = idx
  70. information = []
  71. while count < self.dihedral_type_numbers:
  72. start_idx += 1
  73. if "%FORMAT" in context[start_idx]:
  74. continue
  75. else:
  76. value = list(map(float, context[start_idx].strip().split()))
  77. information.extend(value)
  78. count += len(value)
  79. self.phase_type = information[:self.dihedral_type_numbers]
  80. break
  81. for idx, val in enumerate(context):
  82. if "%FLAG DIHEDRAL_PERIODICITY" in val:
  83. count = 0
  84. start_idx = idx
  85. information = []
  86. while count < self.dihedral_type_numbers:
  87. start_idx += 1
  88. if "%FORMAT" in context[start_idx]:
  89. continue
  90. else:
  91. value = list(map(float, context[start_idx].strip().split()))
  92. information.extend(value)
  93. count += len(value)
  94. self.pn_type = information[:self.dihedral_type_numbers]
  95. break
  96. self.processor(context)
  97. def processor(self, context):
  98. '''processor'''
  99. self.h_atom_a = [0] * self.dihedral_numbers
  100. self.h_atom_b = [0] * self.dihedral_numbers
  101. self.h_atom_c = [0] * self.dihedral_numbers
  102. self.h_atom_d = [0] * self.dihedral_numbers
  103. self.pk = []
  104. self.gamc = []
  105. self.gams = []
  106. self.pn = []
  107. self.ipn = []
  108. for idx, val in enumerate(context):
  109. if "%FLAG DIHEDRALS_INC_HYDROGEN" in val:
  110. count = 0
  111. start_idx = idx
  112. information = []
  113. while count < 5 * self.dihedral_with_hydrogen:
  114. start_idx += 1
  115. if "%FORMAT" in context[start_idx]:
  116. continue
  117. else:
  118. value = list(map(int, context[start_idx].strip().split()))
  119. information.extend(value)
  120. count += len(value)
  121. for i in range(self.dihedral_with_hydrogen):
  122. self.h_atom_a[i] = information[i * 5 + 0] / 3
  123. self.h_atom_b[i] = information[i * 5 + 1] / 3
  124. self.h_atom_c[i] = information[i * 5 + 2] / 3
  125. self.h_atom_d[i] = abs(information[i * 5 + 3] / 3)
  126. tmpi = information[i * 5 + 4] - 1
  127. self.pk.append(self.pk_type[tmpi])
  128. tmpf = self.phase_type[tmpi]
  129. if abs(tmpf - self.CONSTANT_Pi) <= 0.001:
  130. tmpf = self.CONSTANT_Pi
  131. tmpf2 = math.cos(tmpf)
  132. if abs(tmpf2) < 1e-6:
  133. tmpf2 = 0
  134. self.gamc.append(tmpf2 * self.pk[i])
  135. tmpf2 = math.sin(tmpf)
  136. if abs(tmpf2) < 1e-6:
  137. tmpf2 = 0
  138. self.gams.append(tmpf2 * self.pk[i])
  139. self.pn.append(abs(self.pn_type[tmpi]))
  140. self.ipn.append(int(self.pn[i] + 0.001))
  141. break
  142. for idx, val in enumerate(context):
  143. if "%FLAG DIHEDRALS_WITHOUT_HYDROGEN" in val:
  144. count = 0
  145. start_idx = idx
  146. information = []
  147. while count < 5 * (self.dihedral_numbers - self.dihedral_with_hydrogen):
  148. start_idx += 1
  149. if "%FORMAT" in context[start_idx]:
  150. continue
  151. else:
  152. value = list(map(int, context[start_idx].strip().split()))
  153. information.extend(value)
  154. count += len(value)
  155. for i in range(self.dihedral_with_hydrogen, self.dihedral_numbers):
  156. self.h_atom_a[i] = information[(i - self.dihedral_with_hydrogen) * 5 + 0] / 3
  157. self.h_atom_b[i] = information[(i - self.dihedral_with_hydrogen) * 5 + 1] / 3
  158. self.h_atom_c[i] = information[(i - self.dihedral_with_hydrogen) * 5 + 2] / 3
  159. self.h_atom_d[i] = abs(information[(i - self.dihedral_with_hydrogen) * 5 + 3] / 3)
  160. tmpi = information[(i - self.dihedral_with_hydrogen) * 5 + 4] - 1
  161. self.pk.append(self.pk_type[tmpi])
  162. tmpf = self.phase_type[tmpi]
  163. if abs(tmpf - self.CONSTANT_Pi) <= 0.001:
  164. tmpf = self.CONSTANT_Pi
  165. tmpf2 = math.cos(tmpf)
  166. if abs(tmpf2) < 1e-6:
  167. tmpf2 = 0
  168. self.gamc.append(tmpf2 * self.pk[i])
  169. tmpf2 = math.sin(tmpf)
  170. if abs(tmpf2) < 1e-6:
  171. tmpf2 = 0
  172. self.gams.append(tmpf2 * self.pk[i])
  173. self.pn.append(abs(self.pn_type[tmpi]))
  174. self.ipn.append(int(self.pn[i] + 0.001))
  175. break
  176. for i in range(self.dihedral_numbers):
  177. if self.h_atom_c[i] < 0:
  178. self.h_atom_c[i] *= -1