Coverage for gamdpy/configuration/topology.py: 80%

1import numpy as np 

2import numba 

3from numba import cuda 

4 

5class Topology(): 

6 """  

7 This class contains information about the topology, e.g. which bonds, angles and dihedrals are in the system. 

8 """ 

9 

10 def __init__(self, molecule_names=None): 

11 self.bonds = [] 

12 self.angles = [] 

13 self.dihedrals = [] 

14 self.molecules = {} 

15 if molecule_names: 

16 for molecule_name in molecule_names: 

17 self.add_molecule_name(molecule_name) 

18 

19 def add_molecule_name(self, name: str): 

20 self.molecules[name] = [] 

21 

22 # def read() 

23 # def write__to_hdf5() 

24 # ... 

25 

26 def save(self, h5group): 

27 h5group.create_dataset('bonds', data=self.bonds, dtype=np.int32 ) 

28 h5group.create_dataset('angles', data=self.angles, dtype=np.int32 ) 

29 h5group.create_dataset('dihedrals', data=self.dihedrals, dtype=np.int32 ) 

30 

31 h5group.create_group('molecules') 

32 h5group['molecules'].attrs['names'] = list(self.molecules.keys()) # list of names of molecule types 

33 for key in self.molecules.keys(): 

34 h5group['molecules'].create_dataset(key, data=self.molecules[key], dtype=np.int32) 

35 return 

36 

37def bonds_from_positions(positions, cut_off, bond_type): 

38 bonds = [] 

39 for i in range(len(positions)): 

40 for j in range(i): 

41 squared_distance = np.sum( (np.array(positions[j]) - np.array(positions[i]))**2 ) 

42 if squared_distance <= cut_off**2: 

43 bonds.append([j, i, bond_type]) 

44 return bonds 

45 

46def angles_from_bonds(bonds, angle_type): 

47 angles = [] 

48 for bond_index, bond in enumerate(bonds): 

49 for other_bond in bonds[bond_index+1:]: 

50 if bond[0] == other_bond[0]: 

51 angles.append([bond[1], bond[0], other_bond[1], angle_type]) 

52 elif bond[0] == other_bond[1]: 

53 angles.append([bond[1], bond[0], other_bond[0], angle_type]) 

54 elif bond[1] == other_bond[0]: 

55 angles.append([bond[0], bond[1], other_bond[1], angle_type]) 

56 elif bond[1] == other_bond[1]: 

57 angles.append([bond[0], bond[1], other_bond[0], angle_type]) 

58 return angles 

59 

60def dihedrals_from_angles(angles, dihedral_type): 

61 dihedrals = [] 

62 for angle_index, angle in enumerate(angles): 

63 for other_angle in angles[angle_index+1:]: 

64 if angle[1] == other_angle[0] and angle[2] == other_angle[1]: 

65 dihedrals.append([angle[0], angle[1], angle[2], other_angle[2], dihedral_type]) 

66 elif angle[1] == other_angle[2] and angle[2] == other_angle[1]: 

67 dihedrals.append([angle[0], angle[1], angle[2], other_angle[0], dihedral_type]) 

68 elif angle[1] == other_angle[0] and angle[0] == other_angle[1]: 

69 dihedrals.append([angle[2], angle[1], angle[0], other_angle[2], dihedral_type]) 

70 elif angle[1] == other_angle[2] and angle[0] == other_angle[1]: 

71 dihedrals.append([angle[2], angle[1], angle[0], other_angle[0], dihedral_type]) 

72 

73 return dihedrals 

74 

75def molecules_from_bonds(bonds): 

76 molecules = [] 

77 for bond in bonds: 

78 found = False 

79 for molecule in molecules: 

80 if bond[0] in molecule: 

81 found = True 

82 if bond[1] not in molecule: 

83 molecule.append(bond[1]) 

84 break 

85 if bond[1] in molecule: 

86 found = True 

87 if bond[0] not in molecule: 

88 molecule.append(bond[0]) 

89 break 

90 if not found: 

91 molecules.append(bond[0:2]) # If not bonded to existing molecule, the bond is part of a new molecule 

92 return molecules 

93 

94def duplicate_topology(topology, num_molecules): 

95 new_topology = Topology() 

96 assert len(topology.molecules)==1 # Only one type (for now) 

97 for molecule_type in topology.molecules: 

98 assert len(topology.molecules[molecule_type])==1 # Only one molecule 

99 molecule_name = molecule_type 

100 new_topology.add_molecule_name(molecule_name) 

101 particles_per_molecule = len(topology.molecules[molecule_name][0]) 

102 for molecule in range(num_molecules): 

103 first = molecule * particles_per_molecule 

104 for bond in topology.bonds: 

105 new_topology.bonds.append([bond[0] + first, bond[1] + first, bond[2]]) 

106 for angle in topology.angles: 

107 new_topology.angles.append([angle[0] + first, angle[1] + first, angle[2]+ first, angle[3]]) 

108 for dihedral in topology.dihedrals: 

109 new_topology.dihedrals.append(dihedral.copy()) # Copy needed? 

110 for index in range(4): 

111 new_topology.dihedrals[-1][index] += first 

112 new_topology.molecules[molecule_name].append([index + first for index in topology.molecules[molecule_name][0]]) 

113 return new_topology 

114 

115 

116def replicate_topologies(mol_topology_list, num_molecules_each_type_list, mol_types_list, size_molecules_type_list): 

117 """ 

118 

119 Parameters 

120 ---------- 

121 mol_topology_list : list 

122 topology objects, one per moleculke type, to be replicated 

123 num_molecules_each_type_list : list of integers 

124 specifying how many molecules of each type 

125 mol_types_list : list of integers 

126 length should be the total number of molecules, ie the sum of the elements of num_molecules_each_type_list 

127 

128 Returns 

129 ------- 

130 new_topology : topology object representing the full topology of the system of many molecules 

131 

132 """ 

133 new_topology = Topology() 

134 num_molecule_types = len(mol_topology_list) 

135 total_num_molecules = len(mol_types_list) 

136 # check for consistency 

137 assert total_num_molecules == np.sum(np.array(num_molecules_each_type_list)) 

138 

139 for idx in range(num_molecule_types): 

140 assert len(mol_topology_list[idx].molecules) == 1 # only one molecule in each individual topology! 

141 

142 offset = 0 

143 for molecule in range(total_num_molecules): 

144 this_mol_type = mol_types_list[molecule] 

145 num_particles_this_mol = size_molecules_type_list[this_mol_type] 

146 top_this_mol = mol_topology_list[this_mol_type] 

147 for molecule_type in top_this_mol.molecules: 

148 assert len(top_this_mol.molecules[molecule_type])==1 

149 molecule_name = molecule_type 

150 

151 for bond in top_this_mol.bonds: 

152 new_topology.bonds.append([bond[0] + offset, bond[1] + offset, bond[2]]) 

153 for angle in top_this_mol.angles: 

154 new_topology.angles.append([angle[0] + offset, angle[1] + offset, angle[2]+ offset, angle[3]]) 

155 for dihedral in top_this_mol.dihedrals: 

156 new_topology.dihedrals.append(dihedral.copy()) # Copy needed? 

157 for index in range(4): 

158 new_topology.dihedrals[-1][index] += offset 

159 

160 if molecule_name not in new_topology.molecules: 

161 new_topology.add_molecule_name(molecule_name) 

162 new_topology.molecules[molecule_name].append([index + offset for index in top_this_mol.molecules[molecule_name][0]]) 

163 

164 

165 offset += num_particles_this_mol 

166 

167 return new_topology