1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23 from numpy import float64, zeros
24 from warnings import warn
25
26
27 from generic_fns.mol_res_spin import return_molecule, return_residue, return_spin
28 from physical_constants import return_atomic_mass
29 from relax_errors import RelaxError, RelaxNoPdbError
30 from relax_warnings import RelaxWarning
31
32
33
35 """Calculate and return the centre of mass of the structure.
36
37 @keyword atom_id: The molecule, residue, and atom identifier string. Only atoms matching this selection will be used.
38 @type atom_id: str or None
39 @keyword model: Only use a specific model.
40 @type model: int or None
41 @keyword return_mass: A flag which if False will cause only the centre of mass to be returned, but if True will cause the centre of mass and the mass itself to be returned as a tuple.
42 @type return_mass: bool
43 @return: The centre of mass vector, and additionally the mass.
44 @rtype: list of 3 floats (or tuple of a list of 3 floats and one float)
45 """
46
47
48 if not hasattr(cdp, 'structure'):
49 raise RelaxNoPdbError
50
51
52 print("Calculating the centre of mass.")
53
54
55 R = zeros(3, float64)
56
57
58 M = 0.0
59
60
61 for mol_name, res_num, res_name, atom_num, atom_name, element, pos in cdp.structure.atom_loop(atom_id=atom_id, model_num=model, mol_name_flag=True, res_num_flag=True, res_name_flag=True, atom_num_flag=True, atom_name_flag=True, element_flag=True, pos_flag=True):
62
63 id = ''
64
65
66 if mol_name == None:
67 mol_cont = cdp.mol[0]
68 else:
69 id = id + '#' + mol_name
70 mol_cont = return_molecule(id)
71
72
73 if res_name == None and res_num == None:
74 res_cont = mol_cont.res[0]
75 else:
76 id = id + ':' + repr(res_num)
77 res_cont = return_residue(id)
78
79
80 if atom_name == None and atom_num == None:
81 spin_cont = res_cont.spin[0]
82 else:
83 id = id + '@' + repr(atom_num)
84 spin_cont = return_spin(id)
85
86
87 if spin_cont and not spin_cont.select:
88 continue
89
90
91 if element == None:
92 warn(RelaxWarning("Skipping the atom '%s' as the element type cannot be determined." % id))
93 continue
94
95
96 try:
97 mass = return_atomic_mass(element)
98 except RelaxError:
99 warn(RelaxWarning("Skipping the atom '%s' as the element '%s' is unknown." % (id, element)))
100
101
102 M = M + mass
103
104
105 R = R + mass * pos
106
107
108 R = R / M
109
110
111 print(" Total mass: M = " + repr(M))
112 print(" Centre of mass: R = " + repr(R))
113
114
115 if return_mass:
116 return R, M
117 else:
118 return R
119