1 """
2 This modules implements the base class for all the analysis available in nMOLDYN.
3 """
4
5
6 from distutils.version import LooseVersion
7 import getpass
8 import os
9 import platform
10 from random import sample, uniform
11 import re
12 import sys
13 from time import asctime
14 from timeit import default_timer
15
16
17 from Scientific import N
18 from Scientific.Geometry import Tensor, Vector
19
20
21 from MMTK import Atom, AtomCluster, Molecule
22 from MMTK.Collections import Collection
23 from MMTK.NucleicAcids import NucleotideChain
24 from MMTK.ParticleProperties import ParticleScalar
25 from MMTK.PDB import PDBConfiguration
26 from MMTK.Proteins import PeptideChain, Protein
27 from MMTK.Trajectory import Trajectory, TrajectorySet
28 from MMTK import Units
29
30
31 import nMOLDYN
32 from nMOLDYN.Core.Chemistry import belongToAnAmine, belongToAHydroxy, belongToAMethyl, belongToAThiol
33 from nMOLDYN.Core.Error import Error
34 from nMOLDYN.Core.IOFiles import TemporaryFile
35 from nMOLDYN.Core.Logger import LogMessage
36 from nMOLDYN.Core.Mathematics import randomVector
37 from nMOLDYN.Core.Preferences import PREFERENCES
38
39
40 residusChemFamily = {'acidic' : ('Asp','Glu'),
41 'aliphatic' : ('Ile','Leu','Val'),
42 'aromatic' : ('His','Phe','Trp','Tyr'),
43 'basic' : ('Arg','His','Lys'),
44 'charged' : ('Arg','Asp','Glu','His','Lys'),
45 'hydrophobic' : ('Ala','Cys','Cyx','Gly','His','Ile','Leu','Lys','Met','Phe','Thr','Trp','Tyr','Val'),
46 'polar' : ('Arg','Asn','Asp','Cys','Gln','Glu','His','Lys','Ser','Thr','Trp','Tyr'),
47 'small' : ('Ala','Asn','Asp','Cys','Cyx','Gly','Pro','Ser','Thr','Val')}
48
49
50 nmoldyn_package_path = os.path.dirname(os.path.split(__file__)[0])
51
53 """Base class for all analysis defined in nMOLDYN.
54
55 The class Analysis is an abstract-base-class that defines attributes and methods common to all the analysis
56 available in nMOLDYN. To set up an analysis object, use one of its subclass.
57 """
58
59
60 - def __init__(self, parameters = None, statusBar = None):
61 """
62 The constructor.
63
64 @param parameters: a dictionnary that contains parameters of the selected analysis.
65 @type parameters: dict
66
67 @param statusBar: if not None, an instance of nMOLDYN.GUI.Widgets.StatusBar. Will attach a status bar to the
68 selected analysis.
69 @type statusBar: instance of nMOLDYN.GUI.Widgets.StatusBar
70 """
71
72
73 self.chrono = None
74
75 try:
76 self.displayProgressRate = int(PREFERENCES.progress_rate)
77
78 except ValueError:
79 self.displayProgressRate = 10
80
81 if (self.displayProgressRate <= 1) or (self.displayProgressRate >= 100):
82 self.displayProgressRate = 10
83
84
85 if statusBar is None:
86 self.statusBar = None
87 else:
88 self.statusBar = statusBar
89
90 if isinstance(parameters, dict):
91 self.parameters = parameters
92
93 else:
94 self.parameters = None
95
104
105
668
670 """Builds some attributes related to the frame selection string. They will be used to define at which
671 times a given analysis should be run.
672 """
673
674 trajLength = len(self.trajectory)
675
676 if self.first < 0:
677 self.first = 0
678 LogMessage('warning', 'The first step must be an integer >= 1.\n\
679 It will be set to 1 for the running analysis.',['console'])
680
681 if (self.last <= self.first) or (self.last > trajLength):
682 self.last = len(self.trajectory)
683 LogMessage('warning', 'The last step must be an integer in [1,%s].\n\
684 It will be set to %s for the running analysis.' % (trajLength, trajLength), ['console'])
685
686 if (self.skip <= 0) or (self.skip >= trajLength):
687 self.skip = 1
688 LogMessage('warning', 'The skip step must be an integer in [1,%s].\n\
689 It will be set to 1 for the running analysis.' % trajLength,['console'])
690
691 self.frameIndexes = N.arange(self.first, self.last, self.skip)
692 if len(self.frameIndexes) == 0:
693 self.first = 1
694 self.last = trajLength
695 self.skip = 1
696 self.frameIndexes = N.arange(self.first,self.last,self.skip)
697 LogMessage('warning','The frame selection is empty.\
698 It will be set to %s:%s:%s for the running analysis.' % (1,trajLength,1),['console'])
699
700 if hasattr(self.trajectory,'time'):
701 t = self.trajectory.time[self.first:self.last:self.skip]
702
703 else:
704 t = range(trajLength)[self.first:self.last:self.skip]
705
706 if len(t) == 1:
707 self.dt = 1.0
708 else:
709 if t[0] == t[1]:
710 self.dt = 1.0
711 else:
712 self.dt = t[1] - t[0]
713
714 self.times = self.dt * N.arange(len(t))
715
716 self.nFrames = len(self.times)
717
719 """
720 """
721
722 if '_estimate' in self.__class__.__name__:
723 self.preLoadedTraj = None
724 return
725
726 try:
727 LogMessage('info', 'Trying to preload the trajectory. Please wait ...', ['console'])
728 self.preLoadedTraj = {}
729 if structure == 'atom':
730 test = N.zeros((self.nAtoms, self.nFrames, 3), N.Float)
731 del test
732 for atom in self.subset:
733 if differentiation == 0:
734 self.preLoadedTraj[atom] = self.trajectory.readParticleTrajectory(atom, first = self.first, last = self.last, skip = self.skip, variable = 'velocities').array
735 else:
736 self.preLoadedTraj[atom] = self.trajectory.readParticleTrajectory(atom, first = self.first, last = self.last, skip = self.skip).array
737 elif structure == 'frame':
738 test = N.zeros((self.nFrames, self.nAtoms, 3), N.Float)
739 del test
740 for frameIndex in range(self.nFrames):
741 frame = self.frameIndexes[frameIndex]
742 self.preLoadedTraj[frameIndex] = self.trajectory.configuration[frame]
743 LogMessage('info', 'Enough memory to preload the trajectory.', ['console'])
744
745 except:
746 self.preLoadedTraj = None
747 LogMessage('info', 'Not enough memory to preload the trajectory. The analysis will run with runtime %s NetCDF extraction.' % structure, ['console'])
748
750 """Saves the settings of an analysis to an output file.
751
752 @param filename: the name of the output file. If the extension is '.nmi' the output file will be a
753 nMOLDYN input script otherwise the output file will be a nMOLDYN autostart script.
754 @type: filename: string
755 """
756
757 if filename[-4:].lower() == '.nmi':
758
759
760 pScript = open(filename, 'w')
761
762 pScript.write('################################################################\n')
763 pScript.write('# This is an automatically generated python-nMOLDYN run script #\n')
764 pScript.write('################################################################\n\n')
765
766
767 pScript.write(('\nanalysis = %s\n\n') % (self.__class__.__name__))
768
769 pScript.write('# General parameters\n')
770 pScript.write(('version = "%s"\n') % nMOLDYN.__version__)
771 pScript.write(('estimate = "%s"\n\n') % self.parameters['estimate'])
772
773 pScript.write('# Analysis-specific parameters\n')
774
775 for pName in sorted(self.inputParametersNames):
776 pValue = self.parameters[pName]
777
778 if isinstance(pValue, str):
779 pScript.write(('%s = "%s"\n') % (pName, pValue))
780 else:
781 pScript.write(('%s = %s\n') % (pName, pValue))
782
783 pScript.close()
784
785 else:
786
787
788 pScript = open(filename, 'w')
789 pScript.write('#!%s\n' % sys.executable)
790
791
792 pScript.write('################################################################\n')
793 pScript.write('# This is an automatically generated python-nMOLDYN run script #\n')
794 pScript.write('################################################################\n\n')
795
796 pScript.write('from MMTK import *\n')
797 pScript.write(('from %s import %s\n\n') % (self.__module__, self.__class__.__name__))
798
799
800 pScript.write('parameters = {}\n\n')
801
802 pScript.write('# General parameters\n')
803 pScript.write('parameters[\'version\'] = "%s"\n' % nMOLDYN.__version__)
804 pScript.write('parameters[\'estimate\'] = "%s"\n\n' % self.parameters['estimate'])
805
806 pScript.write('# Analysis-specific parameters\n')
807
808 for pName in sorted(self.inputParametersNames):
809 pValue = self.parameters[pName]
810
811 if isinstance(pValue, str):
812 pScript.write(('parameters[\'%s\'] = "%s"\n') % (pName, pValue))
813 else:
814 pScript.write(('parameters[\'%s\'] = %s\n') % (pName, pValue))
815
816
817 pScript.write(('\nanalysis = %s(parameters)') % (self.__class__.__name__))
818
819 pScript.write('\nanalysis.runAnalysis()')
820 pScript.close()
821
822
823 if os.name == 'posix':
824 os.system('chmod u+x %s' % filename)
825
827 """Runs an analysis.
828
829 @return: a dictionnary of the form {'days' : d, 'hours' : h, 'minutes' : m, 'seconds' : s} specifying the
830 time the analysis took in dayx, hours, minutes and seconds.
831 @rtype: dict
832 """
833
834
835 self.jobCounter = 0
836
837
838 self.initialize()
839
840
841 self.buildJobInfo()
842
843
844 if self.version < LooseVersion(vstring = nMOLDYN.__version__):
845
846 LogMessage('warning','The nMOLDYN version number of the script is unknown or deprecated.\n\
847 Some analysis keywords name/values may have changed. Please check in case of troubles.',['file','console'])
848
849
850 self.internalRun()
851
852
853 LogMessage('info','Job finished on '+asctime()+'.\n',['console','file'])
854
855
856 timeTakenForAnalysis = self.analysisTime(self.chrono)
857
858
859
860 if self.estimate:
861
862 LogMessage('info','Estimated time for analysis: %(days)s days %(hours)s hours %(minutes)s minutes %(seconds)s seconds.\n' % timeTakenForAnalysis,['file','console'])
863 else:
864
865 LogMessage('info','Elapsed time for analysis: %(days)s days %(hours)s hours %(minutes)s minutes %(seconds)s seconds.\n' % timeTakenForAnalysis,['file','console'])
866
867 return timeTakenForAnalysis
868
870 """Check the progress of the running analysis and displays periodically on the console and the logfile how far is the analysis.
871 Called each time a step of an analysis loop is achieved.
872
873 @param norm: the maximum number of steps of the analysis.
874 @type: norm: integer
875 """
876
877
878 t = int(100*self.jobCounter/norm)
879
880 oldProgress = (t/self.displayProgressRate)*self.displayProgressRate
881
882
883 self.jobCounter += 1
884
885
886 t = int(100*self.jobCounter/norm)
887 newProgress = (t/self.displayProgressRate)*self.displayProgressRate
888
889 if newProgress != oldProgress:
890 LogMessage('info','Progress rate = %3d %%' % newProgress,['file','console'])
891
892 if self.statusBar is not None:
893 self.statusBar.setValue(t)
894
896 """Display on the console and in the log file the main ifnormation about the analysis to run.
897 """
898
899 self.information = ""
900
901
902 jobCreationTime = asctime()
903
904 self.information += '#'*90 +'\n'
905 self.information += 'Job information for %s analysis.\n' % self.__class__.__name__
906 self.information += '#'*90 +'\n\n'
907
908
909 self.information += 'Job launched on: %s\n\n' % jobCreationTime
910 self.information += 'General informations\n'
911 self.information += '--------------------\n'
912 self.information += 'User: %s\n' % getpass.getuser()
913 self.information += 'OS: %s-%s\n' % (platform.system(), platform.release())
914 self.information += 'Processor: %s\n' % platform.machine()
915 self.information += 'nMOLDYN version: %s\n' % self.parameters['version']
916 self.information += 'Estimate run: %s\n\n' % self.parameters['estimate']
917
918
919 self.information += 'Parameters\n'
920 self.information += '----------\n'
921 self.information += '\n'.join(['%s = %s' % (k,self.parameters[k]) for k in self.parameters if k not in ['version','estimate']])
922 self.information += '\n\n\n'
923
924
925 if hasattr(self, 'subset'):
926 self.information += 'Subset selection\n'
927 self.information += '----------------\n'
928 self.information += 'Number of atoms selected for analysis: %s\n\n' % len(self.subset)
929
930 if hasattr(self, 'deuteration'):
931 self.information += 'Deuteration selection\n'
932 self.information += '---------------------\n'
933 self.information += 'Number of atoms selected for deuteration: %s\n\n' % len(self.deuteration)
934
935 if hasattr(self, 'group'):
936 self.information += 'Group selection\n'
937 self.information += '---------------------\n'
938 self.information += 'Number of groups selected: %s\n\n' % len(self.group)
939
940 if hasattr(self, 'triplet'):
941 self.information += 'Triplet selection\n'
942 self.information += '---------------------\n'
943 self.information += 'Number of triplets selected: %s\n\n' % len(self.triplet)
944
945 if hasattr(self, 'bond'):
946 self.information += 'Bond selection\n'
947 self.information += '---------------------\n'
948 self.information += 'Number of bonds selected: %s\n\n' % len(self.bond)
949
950 self.information += '\n'
951 self.information += 'Job status\n'
952 self.information += '----------\n'
953
954 [LogMessage('info', l, ['file', 'console']) for l in self.information.splitlines()]
955
957 """Converts a time in second in days, hours, minutes and seconds.
958
959 @param time: the time (in seconds) to convert.
960 @type time: integer.
961
962 @return: a dictionnary of the form {'days' : d, 'hours' : h, 'minutes' : m, 'seconds' : s}
963 where d, h, m and s are integers resulting respectively from the conversion of |time| in days, hours, minutes
964 and seconds.
965 @rtype: dict
966 """
967
968 if time is None:
969 return None
970
971 else:
972 converted = {'days' : 0, 'hours' : 0, 'minutes' : 0, 'seconds' : 0}
973
974 converted['days'], time = divmod(time, 86400)
975 converted['hours'], time = divmod(time, 3600)
976 converted['minutes'], converted['seconds'] = divmod(time, 60)
977
978 return converted
979
981 """Returns the weights of |atoms| MMTK collection of |universe| MMTK universe using the
982 weighting scheme |scheme|.
983
984 @param universe: the MMTK universe.
985 @type universe: instance of MMTK.Universe
986
987 @param atoms: the atoms to take into account when defining the weights.
988 @type atoms: instance of MMTK.Collections.Collection
989
990 @param deuter: the hydrogen atoms that will be parametrized as deuterium atoms.
991 @type deuter: instance of MMTK.Collections.Collection
992
993 @param scheme: a string equal to 'equal', 'mass', 'coherent' , 'incoherent' or 'atomicNumber' that specifies the weighting
994 scheme to use.
995 @type scheme: string
996
997 @return: the weights of the selected atoms.
998 @rtype: an instance of MMTK.ParticleProperties.ParticledScalar
999 """
1000
1001 try:
1002 if scheme == 'equal':
1003
1004 weights = ParticleScalar(universe, N.ones(universe.numberOfAtoms(), N.Float))
1005
1006 elif scheme == 'mass':
1007
1008 weights = universe.masses()
1009
1010 for atom in deuter:
1011 if abs(atom._mass - atom.mass_deut) < 1.e-10:
1012 LogMessage('warning','The atom %s is not sensitive to deuteration.' % str(atom),['file'])
1013 else:
1014 weights[atom] = atom.mass_deut
1015
1016 elif scheme == 'atomicNumber':
1017 try:
1018
1019 weights = ParticleScalar(universe)
1020 for at in universe.atomList():
1021 weights[at] = at.atomic_number
1022 except KeyError:
1023 raise Error('One or more element with unknown atomic number in your selection.')
1024
1025 elif scheme == 'coherent':
1026
1027
1028 weights = universe.getAtomScalarArray('b_coherent')
1029
1030
1031 for atom in deuter:
1032 if abs(atom.b_coherent - atom.b_coherent_deut) < 1.e-10:
1033 LogMessage('warning','The atom %s is not sensitive to deuteration.' % str(atom),['file'])
1034 else:
1035 weights[atom] = atom.b_coherent_deut
1036
1037 elif scheme == 'incoherent':
1038
1039
1040 weights = universe.getAtomScalarArray('b_incoherent')
1041
1042
1043 for atom in deuter:
1044 if abs(atom.b_incoherent - atom.b_incoherent_deut) < 1.e-10:
1045 LogMessage('warning','The atom %s is not sensitive to deuteration.' % str(atom),['file'])
1046 else:
1047 weights[atom] = atom.b_incoherent_deut
1048
1049
1050 weights = weights*weights
1051
1052 else:
1053 raise Error('%s weighting scheme is unknown for %s analysis.' % (scheme, self.__class__.__name__))
1054
1055
1056 if atoms != universe:
1057 weights = weights*atoms.booleanMask()
1058
1059
1060 if scheme == 'coherent':
1061 weights = weights/N.sqrt((weights*weights).sumOverParticles())
1062 else:
1063 weights = weights/weights.sumOverParticles()
1064
1065 except:
1066 raise Error('Problem when setting up %s weighting scheme for %s analysis.' % (scheme, self.__class__.__name__))
1067 else:
1068 return weights
1069
1071 """Returns a MMTK collection of atoms that matches |selection| selection string. Used to apply an analysis
1072 to a subset of atoms.
1073
1074 @param universe: the universe on which the selection will be performed.
1075 @type universe: instance of MMTK.Universe
1076
1077 @param selection: the selection string that will define the atoms to select.
1078 @type selection: string
1079
1080 @return: a MMTK Collection of the atoms that matches |selection| selection string.
1081 @rtype: instance of MMTK.Collections.Collection
1082 """
1083
1084 if not isinstance(selection, (Collection, list, tuple, str)):
1085 raise Error('Wrong format for the subset selection. Must be a string or a MMTK Collection' % str(selection))
1086
1087 if isinstance(selection, Collection):
1088 return selection
1089
1090 elif isinstance(selection, (list, tuple)):
1091 try:
1092 subsetSelection = Collection(selection)
1093 return subsetSelection
1094
1095 except:
1096 raise Error('The results of a subset selection must be a MMTK Collection.')
1097
1098
1099 elif selection.lower() == 'all':
1100 return Collection(universe.atomList())
1101
1102 else:
1103
1104
1105 if selection.lower().strip()[0:8] == 'filename':
1106 subsetSelection = self.__parseFileSelectionString(universe, selection, 'subset')
1107
1108 elif selection.lower().strip()[0:10] == 'expression':
1109 subsetSelection = self.__parseExpressionSelectionString(universe, selection, 'subset')
1110
1111
1112 else:
1113 subsetSelection = self.__parseObjectSelectionString(universe, selection)
1114
1115 try:
1116 subsetSelection = Collection(list(subsetSelection))
1117
1118 except:
1119 raise Error('The results of a subset selection must be a MMTK Collection.')
1120
1121 return subsetSelection
1122
1124 """Returns a MMTK collection of atoms that matches |selection| selection string. Used to switch the parameters
1125 of a subset (or all) of hydrogen atoms to the parameters of deuterium in order to simulate deuterated system.
1126
1127 @param universe: the universe on which the selection will be performed.
1128 @type universe: instance of MMTK.Universe
1129
1130 @param selection: the selection string that will define the atoms to select.
1131 @type selection: string
1132
1133 @return: a MMTK Collection of the atoms that matches |selection| selection string.
1134 @rtype: instance of MMTK.Collections.Collection
1135 """
1136
1137 if not isinstance(selection, (Collection, list, tuple, str)):
1138 raise Error('Wrong format for the deuteration selection. Must be a string or a MMTK Collection')
1139
1140 if isinstance(selection, Collection):
1141 return selection
1142
1143 elif isinstance(selection, (list, tuple)):
1144 try:
1145 deuterationSelection = Collection(selection)
1146 return deuterationSelection
1147
1148 except:
1149 raise Error('The results of a deuteration selection must be a MMTK Collection.')
1150
1151 elif selection.lower() == 'no':
1152 return Collection()
1153
1154 elif selection.lower() == 'all':
1155 return Collection([at for at in universe.atomList() if at.type.name == 'hydrogen'])
1156
1157 else:
1158
1159
1160 if selection.lower().strip()[0:8] == 'filename':
1161 deuterationSelection = self.__parseFileSelectionString(universe, selection, 'deuteration')
1162
1163 elif selection.lower().strip()[0:10] == 'expression':
1164 deuterationSelection = self.__parseExpressionSelectionString(universe, selection, 'deuteration')
1165
1166
1167 else:
1168 deuterationSelection = self.__parseObjectSelectionString(universe, selection)
1169
1170 deuterationSelection = [at for at in deuterationSelection if at.type.name == 'hydrogen']
1171
1172 try:
1173 deuterationSelection = Collection(deuterationSelection)
1174
1175 except:
1176 raise Error('The results of a deuteration selection must be a MMTK Collection.')
1177
1178 return deuterationSelection
1179
1181 """Returns a list of MMTK collections where each collection defines a group on which will be applied collectively
1182 an analysis.
1183
1184 @param universe: the universe on which the selection will be performed.
1185 @type universe: instance of MMTK.Universe
1186
1187 @param selection: the selection string that will define the contents of each group.
1188 @type selection: string
1189
1190 @return: a list of MMTK Collection where each collection defines a group..
1191 @rtype: list
1192 """
1193
1194 if not isinstance(selection, (list, tuple, str)):
1195 raise Error('Wrong format for the group selection. It must be a list or tuple of MMTK Collections or a group selection string.')
1196
1197 if isinstance(selection, (list, tuple)):
1198 for el in selection:
1199 if not isinstance(el, Collection):
1200 raise Error('Wrong format for the group selection. It must be a list or tuple of MMTK Collections or a group selection string.')
1201
1202 return list(selection)
1203
1204 elif selection.lower() == 'all':
1205 groupSelection = []
1206 for obj in universe.objectList():
1207 if isinstance(obj,(PeptideChain, Protein, NucleotideChain)):
1208 for r in obj.residues():
1209 groupSelection.append(r.atomList())
1210
1211 else:
1212 groupSelection.append(obj.atomList())
1213
1214 return [Collection(list(g)) for g in groupSelection]
1215
1216 else:
1217
1218 groupSelection = []
1219
1220 if selection.lower().strip()[0:8] == 'filename':
1221 groupSelection = self.__parseFileSelectionString(universe, selection, 'group')
1222
1223 elif selection.lower().strip()[0:10] == 'expression':
1224 groupSelection = self.__parseExpressionSelectionString(universe, selection, 'group')
1225
1226
1227 else:
1228 selectionStringLists = [v.strip() for v in re.split('group\d+', selection) if v != '']
1229
1230 for selString in selectionStringLists:
1231
1232 groupingLevel, sString = [v.strip() for v in re.split('groupinglevel\s+(\w+)\s+', selString) if v != '']
1233
1234 grp = self.__parseObjectSelectionString(universe, sString)
1235
1236 self.__buildGroup(grp, groupSelection, groupingLevel)
1237
1238 try:
1239 groupSelection = [Collection(list(g)) for g in groupSelection]
1240
1241 except:
1242 raise Error('Wrong format for the group selection. It must be a list or tuple of MMTK Collections or a group selection string.')
1243
1244 return groupSelection
1245
1247 """
1248 """
1249
1250 gDict= {}
1251
1252 if gLevel.lower() == 'chain':
1253 for at in g:
1254 chain = at.parent.parent.parent
1255 if gDict.has_key(chain):
1256 gDict[chain].append(at)
1257 else:
1258 gDict[chain] = [at]
1259
1260 elif gLevel.lower() == 'residue':
1261 for at in g:
1262 residue = at.parent.parent
1263 if gDict.has_key(residue):
1264 gDict[residue].append(at)
1265 else:
1266 gDict[residue] = [at]
1267
1268 elif gLevel.lower() in ['molecule', 'cluster', 'nucleicacid', 'protein']:
1269 for at in g:
1270 molecule = at.topLevelChemicalObject()
1271 if gDict.has_key(molecule):
1272 gDict[molecule].append(at)
1273 else:
1274 gDict[molecule] = [at]
1275
1276 elif gLevel.lower() == 'atom':
1277 for at in g:
1278 gDict[at] = at
1279
1280 elif gLevel.lower() == 'amine':
1281 for at in g:
1282 amine = belongToAnAmine(at)
1283 if amine is not None:
1284 if gDict.has_key(amine):
1285 gDict[amine].append(at)
1286 else:
1287 gDict[amine] = [at]
1288
1289 elif gLevel.lower() == 'hydroxy':
1290 for at in g:
1291 hydroxy = belongToAHydroxy(at)
1292 if hydroxy is not None:
1293 if gDict.has_key(hydroxy):
1294 gDict[hydroxy].append(at)
1295 else:
1296 gDict[hydroxy] = [at]
1297
1298 elif gLevel.lower() == 'methyl':
1299 for at in g:
1300 methyl = belongToAMethyl(at)
1301 if methyl is not None:
1302 if gDict.has_key(methyl):
1303 gDict[methyl].append(at)
1304 else:
1305 gDict[methyl] = [at]
1306
1307 elif gLevel.lower() == 'thiol':
1308 for at in g:
1309 thiol = belongToAnHydroxy(at)
1310 if thiol is not None:
1311 if gDict.has_key(thiol):
1312 gDict[thiol].append(at)
1313 else:
1314 gDict[thiol] = [at]
1315
1316 elif gLevel.lower() == 'default':
1317 for at in g:
1318 obj = at.topLevelChemicalObject()
1319
1320 if isinstance(obj, Atom):
1321 gDict[at] = at
1322 elif isinstance(obj, (AtomCluster,Molecule)):
1323 if gDict.has_key(obj):
1324 gDict[obj].append(at)
1325 else:
1326 gDict[obj] = [at]
1327 elif isinstance(obj, (NucleotideChain,PeptideChain,Protein)):
1328 residue = at.parent.parent
1329 if gDict.has_key(residue):
1330 gDict[residue].append(at)
1331 else:
1332 gDict[residue] = [at]
1333
1334 for v in gDict.values():
1335 gSelection.append(set(v))
1336
1338 """Parses a nMOLDYN selection string that starts with 'expression' keyword.
1339
1340 @param universe: the universe on which the selection will be performed.
1341 @type universe: instance of MMTK.Universe
1342
1343 @param selectionString: the selection string to parse.
1344 @type selectionString: string
1345
1346 @param selectionMode: a string being one of 'subset', 'deuteration' or 'group' that will specify which kind of
1347 selection should be performed.
1348 @type selectionMode: string
1349
1350 @return: will depend on |selectionMode|.
1351 @rtype: Python set for |selectionMode| = 'subset' or 'deuteration' and dict for |selectionMode| = 'group'
1352 @return: a set of the atoms that match |selectionString| selection string.
1353 @rtype: set
1354 """
1355
1356 try:
1357 expression = ' '.join(selectionString.split()[1:])
1358 exec(expression)
1359 except:
1360 raise Error('Could not parse the expression of the selection string.')
1361
1362 if not isinstance(selection,list):
1363 raise Error('The expression of the selection string must provide a list.')
1364
1365 if selectionMode in ['subset', 'deuteration']:
1366
1367 selection = set(selection)
1368
1369 elif selectionMode == 'group':
1370 pass
1371
1372 if not selection:
1373 raise Error('The selection associated with\n%s\nselection string gave an empty selection.\nSomething might be wrong with it.' % selectionString)
1374
1375 return selection
1376
1378 """Parses a nMOLDYN Selection file (nms file) in order to perform a subset, a deuteration or a group
1379 selection.
1380
1381 @param universe: the universe on which the selection will be performed.
1382 @type universe: instance of MMTK.Universe
1383
1384 @param selectionString: the selection string to parse.
1385 @type selectionString: string
1386
1387 @param selectionMode: a string being one of 'subset', 'deuteration' or 'group' that will specify which kind of
1388 selection should be performed.
1389 @type selectionMode: string
1390
1391 @return: will depend on |selectionMode|.
1392 @rtype: Python set for |selectionMode| = 'subset' or 'deuteration' and dict for |selectionMode| = 'group'
1393 """
1394
1395
1396 try:
1397 filename = selectionString.split()[1].strip()
1398 nmsFile = open(filename, 'r')
1399 exec nmsFile
1400 nmsFile.close()
1401 except:
1402 raise Error('Unable to open the file %s.' % filename)
1403
1404 try:
1405 pdb = pdb.strip()
1406 if not os.path.exists(pdb):
1407 raise Error('The pdb file %s does not exist.' % pdb)
1408
1409 except NameError:
1410 raise Error('The "pdb" field was not set in the %s nms file.' % filename)
1411
1412
1413 if selectionMode in ['subset', 'deuteration']:
1414 try:
1415 selectedAtoms = [int(v) for v in eval(selectionMode)]
1416 except:
1417 raise Error('The %s selection field must be of the form %s = [int,int,int...]' % (selectionMode,selectionMode))
1418
1419 selection = set()
1420
1421 elif selectionMode == 'group':
1422 try:
1423 selectedAtoms = [[int(vv) for vv in v] for v in eval(selectionMode)]
1424
1425 except :
1426 raise Error('The group selection field must be of the form group = [[int,int,...],[int,int,...],[int,int,...]...]')
1427
1428 selection = [set()]*len(selectedAtoms)
1429
1430 else:
1431 raise Error('Unknown selection type: %s.' % selectionMode)
1432
1433 try:
1434 atomPos = []
1435
1436 for atom in universe.atomList():
1437 atomPos.append([atom, [round(v,2) for v in self.trajectory.configuration[0][atom]]])
1438
1439
1440 pdb = PDBConfiguration(pdb)
1441
1442 if selectionMode in ['subset','deuteration']:
1443
1444
1445 for res in pdb.residues:
1446
1447 for at in res.atom_list:
1448
1449 if at.properties['serial_number'] in selectedAtoms:
1450 pos = [round(v,2) for v in at.position]
1451
1452 for atom in atomPos:
1453
1454 if atom[1] == pos:
1455
1456 selection.add(atom[0])
1457 break
1458
1459 elif selectionMode == 'group':
1460
1461
1462 for res in pdb.residues:
1463
1464 for at in res.atom_list:
1465
1466 for p in selectedAtoms:
1467
1468 if at.properties['serial_number'] in p:
1469 pos = [round(v,2) for v in at.position]
1470
1471 for atom in atomPos:
1472
1473 if atom[1] == pos:
1474 selection[selectedAtoms.index(p)].add(atom[0])
1475 break
1476
1477 except:
1478 raise Error('Could not parse properly for selection the PDB file associated with %s selection file.' % filename)
1479
1480 if not selection:
1481 raise Error('The selection from file %s gave an empty selection. Something might be wrong with it.' % filename)
1482
1483 return selection
1484
1486 """Parses a selection string.
1487
1488 @param universe: the universe on which the selection will be performed.
1489 @type universe: instance of MMTK.Universe
1490
1491 @param selectionString: the selection strin to parse.
1492 @type selectionString: string
1493
1494 @return: a set of the atoms that match |selectionString| selection string.
1495 @rtype: set
1496 """
1497
1498 selStr = re.sub('\(',' ( ',selectionString)
1499 selStr = re.sub('\)',' ) ',selStr)
1500
1501 l = [v.strip() for v in re.split('\s+|,',selStr) if v.strip()]
1502
1503 processedList = []
1504 selectionValue = []
1505
1506 try:
1507 while True:
1508
1509 if not l:
1510 break
1511
1512 l0 = l[0].lower()
1513 if l0 == 'objectname':
1514 objectName = l[1]
1515 del l[0:2]
1516 objectClass = universe.nmoldyncontents[objectName]['objectclass'].lower()
1517 continue
1518
1519 elif l0 in universe.nmoldyncontents[objectName].keys():
1520 selectionKeyword = l0
1521 del l[0]
1522
1523 while True:
1524
1525 ll0 = l[0].lower()
1526 if ll0 == 'objectname':
1527 raise
1528
1529 elif ll0 in universe.nmoldyncontents[objectName].keys():
1530 raise
1531
1532 elif ll0 in ['and','or',')']:
1533 if not selectionValue:
1534 raise
1535 indexes = self.__retrieveAtomIndexes(universe, objectClass, objectName, selectionKeyword, selectionValue)
1536 processedList.append(str(set(indexes)))
1537 selectionValue = []
1538 break
1539
1540 elif ll0 == '(':
1541 raise
1542
1543 else:
1544 selectionValue.append(ll0)
1545 del l[0]
1546 if not l:
1547 indexes = self.__retrieveAtomIndexes(universe, objectClass, objectName, selectionKeyword, selectionValue)
1548 processedList.append(str(set(indexes)))
1549 selectionValue = []
1550 break
1551
1552 elif l0 in ['and','or','(']:
1553 processedList.append(l0)
1554 del l[0]
1555 if not l:
1556 raise
1557
1558 elif l0 == ')':
1559 processedList.append(l0)
1560 del l[0]
1561
1562 else:
1563 raise
1564
1565 except:
1566 raise Error('Badly formatted selection string.')
1567
1568 processedString = ''.join(processedList)
1569 processedString = processedString.replace('and', '&')
1570 processedString = processedString.replace('or', '|')
1571
1572
1573 selection = set([universe.indexToAtoms[ind] for ind in eval(processedString)])
1574
1575 if not selection:
1576 raise Error('The selection associated with\n%s\nselection string gave an empty selection.\nSomething might be wrong with it.' % selectionString)
1577
1578 return selection
1579
1581 """Retrieves the MMTK index of the atoms matching |objectClass| MMTK chemical object class,
1582 |objectName| nMOLDYN name, |selectionKeyword| selection keyword and |selectionValue| selection value.
1583
1584 @param universe: the universe on which the selection will be performed.
1585 @type universe: instance of MMTK.Universe
1586
1587 @param objectClass: the MMTK chemical object to match.
1588 @type objectClass: string
1589
1590 @param objectName: the nMOLDYN name to match.
1591 @type objectName: string
1592
1593 @param selectionKeyword: the selection keyword to match.
1594 @type selectionKeyword: string
1595
1596 @param selectionValue: the selection value to match.
1597 @type selectionValue: string
1598
1599 @return: a list of MMTK indexes of the selected atoms.
1600 @rtype: list
1601 """
1602
1603 if objectClass == 'allclass':
1604
1605 indexes = self.__allClassParser(universe, selectionKeyword, selectionValue)
1606
1607 elif objectClass == 'atom':
1608
1609 indexes = self.__atomParser(universe, objectName.lower(), selectionKeyword, selectionValue)
1610
1611 elif objectClass == 'atomcluster':
1612
1613 indexes = self.__atomClusterParser(universe, objectName.lower(), selectionKeyword, selectionValue)
1614
1615 elif objectClass == 'molecule':
1616
1617 indexes = self.__moleculeParser(universe, objectName.lower(), selectionKeyword, selectionValue)
1618
1619 elif objectClass == 'nucleotidechain':
1620
1621 indexes = self.__nucleotideChainParser(universe, objectName.lower(), selectionKeyword, selectionValue)
1622
1623 elif objectClass == 'peptidechain':
1624
1625 indexes = self.__peptideChainParser(universe, objectName.lower(), selectionKeyword, selectionValue)
1626
1627 elif objectClass == 'protein':
1628
1629 indexes = self.__proteinParser(universe, objectName.lower(), selectionKeyword, selectionValue)
1630
1631 else:
1632 raise Error('The chemical class %s is not defined in MMTK.' % objectClass)
1633
1634 return indexes
1635
1637 """Retrieves the MMTK indexes of the atoms whose nMOLDYN name is '*' matching |selectionKeyword|
1638 selection keyword and |selectionValue| selection value.
1639
1640 @param universe: the universe on which the selection will be performed.
1641 @type universe: instance of MMTK.Universe
1642
1643 @param selectionKeyword: the selection keyword to match.
1644 @type selectionKeyword: string
1645
1646 @param selectionValue: the selection value to match.
1647 @type selectionValue: string
1648
1649 @return: a list of MMTK indexes of the selected atoms.
1650 @rtype: list
1651 """
1652
1653 selectedObjects = universe.objectList()
1654
1655 selection = []
1656
1657 if selectionKeyword == 'atomelement':
1658 if '*' in selectionValue:
1659 for obj in selectedObjects:
1660 selection.extend([at.index for at in obj.atomList()])
1661 else:
1662 selection = []
1663 for v in selectionValue:
1664 for obj in selectedObjects:
1665 selection.extend([at.index for at in obj.atomList() if at.type.name.strip().lower() == v])
1666
1667 return selection
1668
1669 - def __atomParser(self, universe, objectName, selectionKeyword, selectionValue):
1670 """Retrieves the MMTK index of the atoms whose MMTK chemical object class is 'Atom' matching
1671 |objectName| nMOLDYN name, |selectionKeyword| selection keyword and |selectionValue| selection value.
1672
1673 @param universe: the universe on which the selection will be performed.
1674 @type universe: instance of MMTK.Universe
1675
1676 @param objectName: the nMOLDYN name to match.
1677 @type objectName: string
1678
1679 @param selectionKeyword: the selection keyword to match.
1680 @type selectionKeyword: string
1681
1682 @param selectionValue: the selection value to match.
1683 @type selectionValue: string
1684
1685 @return: a list of MMTK indexes of the selected atoms.
1686 @rtype: list
1687 """
1688
1689 selection = []
1690
1691 if selectionKeyword == 'name':
1692
1693
1694 if '*' in selectionValue:
1695 selection = [obj.index for obj in universe.objectList() if obj.nmoldynname.strip().lower() == objectName]
1696
1697
1698
1699 else:
1700 selection = []
1701 for v in selectionValue:
1702 selection.extend([obj.index for obj in universe.objectList() if obj.name.strip().lower() == v])
1703
1704 return selection
1705
1707 """Retrieves the MMTK index of the atoms whose MMTK chemical object class is 'AtomCluster' matching
1708 |objectName| nMOLDYN name, |selectionKeyword| selection keyword and |selectionValue| selection value.
1709
1710 @param universe: the universe on which the selection will be performed.
1711 @type universe: instance of MMTK.Universe
1712
1713 @param objectName: the nMOLDYN name to match.
1714 @type objectName: string
1715
1716 @param selectionKeyword: the selection keyword to match.
1717 @type selectionKeyword: string
1718
1719 @param selectionValue: the selection value to match.
1720 @type selectionValue: string
1721
1722 @return: a list of MMTK indexes of the selected atoms.
1723 @rtype: list
1724 """
1725
1726 selection = []
1727
1728
1729
1730 if selectionKeyword == 'name':
1731
1732
1733 if '*' in selectionValue:
1734 selectedObjects = set([obj for obj in universe.objectList() if obj.nmoldynname.strip().lower() == objectName])
1735
1736
1737
1738 else:
1739 selectedObjects = set()
1740 for v in selectionValue:
1741 selectedObjects.update([obj for obj in universe.objectList() if obj.name.strip().lower() == v])
1742
1743 for obj in selectedObjects:
1744 selection.extend([at.index for at in obj.atomList()])
1745
1746 else:
1747 selectedObjects = set([obj for obj in universe.objectList() if obj.nmoldynname.strip().lower() == objectName])
1748
1749 if selectionKeyword == 'atomname':
1750 for v in selectionValue:
1751 if v == '*':
1752 for obj in selectedObjects:
1753 selection.extend([at.index for at in obj.atomList()])
1754 else:
1755 for obj in selectedObjects:
1756 selection.extend([at.index for at in obj.atomList() if at.fullName().strip().lower() == v])
1757
1758 elif selectionKeyword == 'atomelement':
1759 for v in selectionValue:
1760 if v == '*':
1761 for obj in selectedObjects:
1762 selection.extend([at.index for at in obj.atomList()])
1763 else:
1764 for obj in selectedObjects:
1765 selection.extend([at.index for at in obj.atomList() if at.type.name.strip().lower() == v])
1766
1767 return selection
1768
1769 - def __moleculeParser(self, universe, objectName, selectionKeyword, selectionValue):
1770 """Retrieves the MMTK index of the atoms whose MMTK chemical object class is 'Molecule' matching
1771 |objectName| nMOLDYN name, |selectionKeyword| selection keyword and |selectionValue| selection value.
1772
1773 @param universe: the universe on which the selection will be performed.
1774 @type universe: instance of MMTK.Universe
1775
1776 @param objectName: the nMOLDYN name to match.
1777 @type objectName: string
1778
1779 @param selectionKeyword: the selection keyword to match.
1780 @type selectionKeyword: string
1781
1782 @param selectionValue: the selection value to match.
1783 @type selectionValue: string
1784
1785 @return: a list of MMTK indexes of the selected atoms.
1786 @rtype: list
1787 """
1788
1789 selection = []
1790
1791
1792 if selectionKeyword == 'name':
1793
1794
1795 if '*' in selectionValue:
1796 selectedObjects = set([obj for obj in universe.objectList() if obj.nmoldynname.strip().lower() == objectName])
1797
1798
1799
1800 else:
1801 selectedObjects = set()
1802 for v in selectionValue:
1803 selectedObjects.update([obj for obj in universe.objectList() if obj.name.strip().lower() == v])
1804 for obj in selectedObjects:
1805 selection.extend([at.index for at in obj.atomList()])
1806
1807 else:
1808 selectedObjects = set([obj for obj in universe.objectList() if obj.nmoldynname.strip().lower() == objectName])
1809
1810 if selectionKeyword == 'atomname':
1811 for v in selectionValue:
1812 if v == '*':
1813 for obj in selectedObjects:
1814 selection.extend([at.index for at in obj.atomList()])
1815 else:
1816 for obj in selectedObjects:
1817 selection.extend([at.index for at in obj.atomList() if at.fullName().strip().lower() == v])
1818
1819 elif selectionKeyword == 'atomelement':
1820 for v in selectionValue:
1821 if v == '*':
1822 for obj in selectedObjects:
1823 selection.extend([at.index for at in obj.atomList()])
1824 else:
1825 for obj in selectedObjects:
1826 selection.extend([at.index for at in obj.atomList() if at.type.name.strip().lower() == v])
1827
1828 elif selectionKeyword == 'chemfragment':
1829 for v in selectionValue:
1830
1831 if v == 'amine':
1832 for obj in selectedObjects:
1833 nitrogens = [at for at in obj.atomList() if at.type.name.strip().lower() == 'nitrogen']
1834 for n in nitrogens:
1835 neighbours = n.bondedTo()
1836 hydrogens = [neigh for neigh in neighbours if neigh.type.name.strip().lower() == 'hydrogen']
1837
1838 if len(hydrogens) == 2:
1839 selection.extend([n.index] + [h.index for h in hydrogens])
1840
1841 elif (len(hydrogens) == 3) and (obj.numberOfAtoms() == 4):
1842 selection.extend([n.index] + [h.index for h in hydrogens])
1843
1844 elif v == 'hydroxy':
1845 for obj in selectedObjects:
1846 oxygens = [at for at in obj.atomList() if at.type.name.strip().lower() == 'oxygen']
1847 for o in oxygens:
1848 neighbours = o.bondedTo()
1849 hydrogens = [neigh for neigh in neighbours if neigh.type.name.strip().lower() == 'hydrogen']
1850
1851 if len(hydrogens) == 1:
1852 selection.extend([o.index] + [h.index for h in hydrogens])
1853
1854 elif (len(hydrogens) == 2) and (obj.numberOfAtoms() == 3):
1855 selection.extend([o.index] + [h.index for h in hydrogens])
1856
1857 elif v == 'methyl':
1858 for obj in selectedObjects:
1859 carbons = [at for at in obj.atomList() if at.type.name.strip().lower() == 'carbon']
1860 for c in carbons:
1861 neighbours = c.bondedTo()
1862 hydrogens = [neigh for neigh in neighbours if neigh.type.name.strip().lower() == 'hydrogen']
1863
1864 if len(hydrogens) == 3:
1865 selection.extend([c.index] + [h.index for h in hydrogens])
1866
1867 elif (len(hydrogens) == 4) and (obj.numberOfAtoms() == 5):
1868 selection.extend([c.index] + [h.index for h in hydrogens])
1869
1870 elif v == 'thiol':
1871 for obj in selectedObjects:
1872 sulphurs = [at for at in obj.atomList() if at.type.name.strip().lower() in ['sulphur','sulfur']]
1873 for s in sulphurs:
1874 neighbours = s.bondedTo()
1875 hydrogens = [neigh for neigh in neighbours if neigh.type.name.strip().lower() == 'hydrogen']
1876
1877 if len(hydrogens) == 1:
1878 selection.extend([o.index] + [h.index for h in hydrogens])
1879
1880 elif (len(hydrogens) == 2) and (obj.numberOfAtoms() == 3):
1881 selection.extend([s.index] + [h.index for h in hydrogens])
1882
1883 return selection
1884
1886 """Retrieves the MMTK index of the atoms whose MMTK chemical object class is 'NucleotideChain' matching
1887 |objectName| nMOLDYN name, |selectionKeyword| selection keyword and |selectionValue| selection value.
1888
1889 @param universe: the universe on which the selection will be performed.
1890 @type universe: instance of MMTK.Universe
1891
1892 @param objectName: the nMOLDYN name to match.
1893 @type objectName: string
1894
1895 @param selectionKeyword: the selection keyword to match.
1896 @type selectionKeyword: string
1897
1898 @param selectionValue: the selection value to match.
1899 @type selectionValue: string
1900
1901 @return: a list of MMTK indexes of the selected atoms.
1902 @rtype: list
1903 """
1904
1905 selection = []
1906
1907
1908
1909 if selectionKeyword == 'name':
1910
1911
1912 if '*' in selectionValue:
1913 selectedObjects = set([obj for obj in universe.objectList() if obj.nmoldynname.strip().lower() == objectName])
1914
1915
1916
1917 else:
1918 selectedObjects = set([])
1919 for v in selectionValue:
1920 selectedObjects.update([obj for obj in universe.objectList() if obj.name.strip().lower() == v])
1921
1922 for obj in selectedObjects:
1923 selection.extend([at.index for at in obj.atomList()])
1924
1925 else:
1926 selectedObjects = set([obj for obj in universe.objectList() if obj.nmoldynname.strip().lower() == objectName])
1927
1928 if selectionKeyword == 'atomname':
1929 for v in selectionValue:
1930 if v == '*':
1931 for obj in selectedObjects:
1932 selection.extend([at.index for at in obj.atomList()])
1933 else:
1934 for obj in selectedObjects:
1935 selection.extend([at.index for at in obj.atomList() if at.fullName().strip().lower() == v])
1936
1937 elif selectionKeyword == 'atomtype':
1938 for v in selectionValue:
1939 if v == '*':
1940 for obj in selectedObjects:
1941 selection.extend([at.index for at in obj.atomList()])
1942
1943 else:
1944 for obj in selectedObjects:
1945 selection.extend([at.index for at in obj.atomList() if at.name.strip().lower() == v])
1946
1947 elif selectionKeyword == 'atomelement':
1948 for v in selectionValue:
1949 if v == '*':
1950 for obj in selectedObjects:
1951 selection.extend([at.index for at in obj.atomList()])
1952 else:
1953 for obj in selectedObjects:
1954 selection.extend([at.index for at in obj.atomList() if at.type.name.strip().lower() == v])
1955
1956 elif selectionKeyword == 'nuclname':
1957 for v in selectionValue:
1958 if v == '*':
1959 for obj in selectedObjects:
1960 selection.extend([at.index for at in obj.atomList()])
1961 else:
1962 for obj in selectedObjects:
1963 for nucl in obj.residues():
1964 if nucl.fullName().strip().lower() == v:
1965 selection.extend([at.index for at in nucl.atomList()])
1966
1967 elif selectionKeyword == 'nucltype':
1968 for v in selectionValue:
1969 if v == '*':
1970 for obj in selectedObjects:
1971 selection.extend([at.index for at in obj.atomList()])
1972 else:
1973 for obj in selectedObjects:
1974 for nucl in obj.residues():
1975 if nucl.symbol.strip().lower() == v:
1976 selection.extend([at.index for at in nucl.atomList()])
1977
1978 elif selectionKeyword == 'misc':
1979 for v in selectionValue:
1980 if v == 'bases':
1981 for obj in selectedObjects:
1982 selection.extend([at.index for at in obj.bases().atomList()])
1983
1984 elif v == 'backbone':
1985 for obj in selectedObjects:
1986 selection.extend([at.index for at in obj.backbone().atomList()])
1987
1988 return selection
1989
1991 """Retrieves the MMTK index of the atoms whose MMTK chemical object class is 'PeptideChain' matching
1992 |objectName| nMOLDYN name, |selectionKeyword| selection keyword and |selectionValue| selection value.
1993
1994 @param universe: the universe on which the selection will be performed.
1995 @type universe: instance of MMTK.Universe
1996
1997 @param objectName: the nMOLDYN name to match.
1998 @type objectName: string
1999
2000 @param selectionKeyword: the selection keyword to match.
2001 @type selectionKeyword: string
2002
2003 @param selectionValue: the selection value to match.
2004 @type selectionValue: string
2005
2006 @return: a list of MMTK indexes of the selected atoms.
2007 @rtype: list
2008 """
2009
2010 selection = []
2011
2012
2013 if selectionKeyword == 'name':
2014
2015
2016
2017 if '*' in selectionValue:
2018 selectedObjects = set([obj for obj in universe.objectList() if obj.nmoldynname.strip().lower() == objectName])
2019
2020
2021
2022 else:
2023 selectedObjects = set([])
2024 for v in selectionValue:
2025 selectedObjects.update([obj.index for obj in universe.objectList() if obj.name.strip().lower() == v])
2026
2027 for obj in selectedObjects:
2028 selection.extend([at.index for at in obj.atomList()])
2029 else:
2030 selectedObjects = set([obj for obj in universe.objectList() if obj.nmoldynname.strip().lower() == objectName])
2031
2032 if selectionKeyword == 'atomname':
2033 for v in selectionValue:
2034 if v == '*':
2035 for obj in selectedObjects:
2036 selection.extend([at.index for at in obj.atomList()])
2037 else:
2038 for obj in selectedObjects:
2039 selection.extend([at.index for at in obj.atomList() if at.fullName().strip().lower() == v])
2040
2041 elif selectionKeyword == 'atomtype':
2042 for v in selectionValue:
2043 if v == '*':
2044 for obj in selectedObjects:
2045 selection.extend([at.index for at in obj.atomList()])
2046 else:
2047 for obj in selectedObjects:
2048 selection.extend([at.index for at in obj.atomList() if at.name.strip().lower() == v])
2049
2050 elif selectionKeyword == 'atomelement':
2051 for v in selectionValue:
2052 if v == '*':
2053 for obj in selectedObjects:
2054 selection.extend([at.index for at in obj.atomList()])
2055 else:
2056 for obj in selectedObjects:
2057 selection.extend([at.index for at in obj.atomList() if at.type.name.strip().lower() == v])
2058
2059 elif selectionKeyword == 'chemfragment':
2060 for v in selectionValue:
2061 if v == 'amine':
2062 for obj in selectedObjects:
2063 nitrogens = [at for at in obj.atomList() if at.type.name.strip().lower() == 'nitrogen']
2064 for n in nitrogens:
2065 neighbours = n.bondedTo()
2066 hydrogens = [neigh for neigh in neighbours if neigh.type.name.strip().lower() == 'hydrogen']
2067
2068 if len(hydrogens) == 2:
2069 selection.extend([n.index] + [h.index for h in hydrogens])
2070
2071 elif (len(hydrogens) == 3) and (obj.numberOfAtoms() == 4):
2072 selection.extend([n.index] + [h.index for h in hydrogens])
2073
2074 elif v == 'c_alphas':
2075 for obj in selectedObjects:
2076 selection.extend([at.index for at in obj.atomList() if at.name.strip().lower() == 'c_alpha'])
2077
2078 elif v == 'hydroxy':
2079 for obj in selectedObjects:
2080 oxygens = [at for at in obj.atomList() if at.type.name.strip().lower() == 'oxygen']
2081 for o in oxygens:
2082 neighbours = o.bondedTo()
2083 hydrogens = [neigh for neigh in neighbours if neigh.type.name.strip().lower() == 'hydrogen']
2084
2085 if len(hydrogens) == 1:
2086 selection.extend([o.index] + [h.index for h in hydrogens])
2087
2088 elif (len(hydrogens) == 2) and (obj.numberOfAtoms() == 3):
2089 selection.extend([o.index] + [h.index for h in hydrogens])
2090
2091 elif v == 'methyl':
2092 for obj in selectedObjects:
2093 carbons = [at for at in obj.atomList() if at.type.name.strip().lower() == 'carbon']
2094 for c in carbons:
2095 neighbours = c.bondedTo()
2096 hydrogens = [neigh for neigh in neighbours if neigh.type.name.strip().lower() == 'hydrogen']
2097
2098 if len(hydrogens) == 3:
2099 selection.extend([c.index] + [h.index for h in hydrogens])
2100
2101 elif (len(hydrogens) == 4) and (obj.numberOfAtoms() == 5):
2102 selection.extend([c.index] + [h.index for h in hydrogens])
2103
2104 elif v == 'thiol':
2105 for obj in selectedObjects:
2106 sulphurs = [at for at in obj.atomList() if at.type.name.strip().lower() in ['sulphur', 'sulfur']]
2107 for s in sulphurs:
2108 neighbours = s.bondedTo()
2109 hydrogens = [neigh for neigh in neighbours if neigh.type.name.strip().lower() == 'hydrogen']
2110
2111 if len(hydrogens) == 1:
2112 selection.extend([o.index] + [h.index for h in hydrogens])
2113
2114 elif (len(hydrogens) == 2) and (obj.numberOfAtoms() == 3):
2115 selection.extend([s.index] + [h.index for h in hydrogens])
2116
2117 elif selectionKeyword == 'resname':
2118 for v in selectionValue:
2119 if v == '*':
2120 for obj in selectedObjects:
2121 selection.extend([at.index for at in obj.atomList()])
2122 else:
2123 for obj in selectedObjects:
2124 for res in obj.residues():
2125 if res.fullName().strip().lower() == v:
2126 selection.extend([at.index for at in res.atomList()])
2127
2128 elif selectionKeyword == 'restype':
2129 for v in selectionValue:
2130 if v == '*':
2131 for obj in selectedObjects:
2132 selection.extend([at.index for at in obj.atomList()])
2133 else:
2134 for obj in selectedObjects:
2135 for res in obj.residues():
2136 if res.symbol.strip().lower() == v:
2137 selection.extend([at.index for at in res.atomList()])
2138
2139 elif selectionKeyword == 'resclass':
2140 for v in selectionValue:
2141 for obj in selectedObjects:
2142 for res in obj.residues():
2143 if res.symbol.strip().lower() in residusChemFamily[v]:
2144 selection.extend([at.index for at in res.atomList()])
2145
2146 elif selectionKeyword == 'misc':
2147 for v in selectionValue:
2148 if v == 'sidechains':
2149 for obj in selectedObjects:
2150 selection.extend([at.index for at in obj.sidechains().atomList()])
2151
2152 elif v == 'backbone':
2153 for obj in selectedObjects:
2154 selection.extend([at.index for at in obj.backbone().atomList()])
2155
2156 return selection
2157
2158 - def __proteinParser(self, universe, objectName, selectionKeyword, selectionValue):
2159 """Retrieves the MMTK index of the atoms whose MMTK chemical object class is 'Protein' matching
2160 |objectName| nMOLDYN name, |selectionKeyword| selection keyword and |selectionValue| selection value.
2161
2162 @param universe: the universe on which the selection will be performed.
2163 @type universe: instance of MMTK.Universe
2164
2165 @param objectName: the nMOLDYN name to match.
2166 @type objectName: string
2167
2168 @param selectionKeyword: the selection keyword to match.
2169 @type selectionKeyword: string
2170
2171 @param selectionValue: the selection value to match.
2172 @type selectionValue: string
2173
2174 @return: a list of MMTK indexes of the selected atoms.
2175 @rtype: list
2176 """
2177
2178 selection = []
2179
2180
2181 if selectionKeyword == 'name':
2182
2183
2184 if '*' in selectionValue:
2185 selectedObjects = set([obj for obj in universe.objectList() if obj.nmoldynname.strip().lower() == objectName])
2186
2187
2188 else:
2189 selectedObjects = set([])
2190 for v in selectionValue:
2191 selectedObjects.update([obj.index for obj in universe.objectList() if obj.name.strip().lower() == v])
2192
2193 for obj in selectedObjects:
2194 selection.extend([at.index for at in obj.atomList()])
2195 else:
2196 selectedObjects = set([obj for obj in universe.objectList() if obj.nmoldynname.strip().lower() == objectName])
2197
2198
2199 if selectionKeyword == 'atomname':
2200 for v in selectionValue:
2201 if v == '*':
2202 for obj in selectedObjects:
2203 selection.extend([at.index for at in obj.atomList()])
2204 else:
2205 for obj in selectedObjects:
2206 selection.extend([at.index for at in obj.atomList() if at.fullName().strip().lower() == v])
2207
2208 elif selectionKeyword == 'atomtype':
2209 for v in selectionValue:
2210 if v == '*':
2211 for obj in selectedObjects:
2212 selection.extend([at.index for at in obj.atomList()])
2213 else:
2214 for obj in selectedObjects:
2215 selection.extend([at.index for at in obj.atomList() if at.name.strip().lower() == v])
2216
2217 elif selectionKeyword == 'atomelement':
2218 for v in selectionValue:
2219 if v == '*':
2220 for obj in selectedObjects:
2221 selection.extend([at.index for at in obj.atomList()])
2222 else:
2223 for obj in selectedObjects:
2224 selection.extend([at.index for at in obj.atomList() if at.type.name.strip().lower() == v])
2225
2226 elif selectionKeyword == 'chemfragment':
2227 for v in selectionValue:
2228
2229 if v == 'amine':
2230 for obj in selectedObjects:
2231 nitrogens = [at for at in obj.atomList() if at.type.name.strip().lower() == 'nitrogen']
2232 for n in nitrogens:
2233 neighbours = n.bondedTo()
2234 hydrogens = [neigh for neigh in neighbours if neigh.type.name.strip().lower() == 'hydrogen']
2235
2236 if len(hydrogens) == 2:
2237 selection.extend([n.index] + [h.index for h in hydrogens])
2238
2239 elif (len(hydrogens) == 3) and (obj.numberOfAtoms() == 4):
2240 selection.extend([n.index] + [h.index for h in hydrogens])
2241
2242 elif v == 'c_alphas':
2243 for obj in selectedObjects:
2244 selection.extend([at.index for at in obj.atomList() if at.name.strip().lower() == 'c_alpha'])
2245
2246 elif v == 'hydroxy':
2247 for obj in selectedObjects:
2248 oxygens = [at for at in obj.atomList() if at.type.name.strip().lower() == 'oxygen']
2249 for o in oxygens:
2250 neighbours = o.bondedTo()
2251 hydrogens = [neigh for neigh in neighbours if neigh.type.name.strip().lower() == 'hydrogen']
2252
2253 if len(hydrogens) == 1:
2254 selection.extend([o.index] + [h.index for h in hydrogens])
2255
2256 elif (len(hydrogens) == 2) and (obj.numberOfAtoms() == 3):
2257 selection.extend([o.index] + [h.index for h in hydrogens])
2258
2259 elif v == 'methyl':
2260 for obj in selectedObjects:
2261 carbons = [at for at in obj.atomList() if at.type.name.strip().lower() == 'carbon']
2262 for c in carbons:
2263 neighbours = c.bondedTo()
2264 hydrogens = [neigh for neigh in neighbours if neigh.type.name.strip().lower() == 'hydrogen']
2265
2266 if len(hydrogens) == 3:
2267 selection.extend([c.index] + [h.index for h in hydrogens])
2268
2269 elif (len(hydrogens) == 4) and (obj.numberOfAtoms() == 5):
2270 selection.extend([c.index] + [h.index for h in hydrogens])
2271
2272 elif v == 'thiol':
2273 for obj in selectedObjects:
2274 sulphurs = [at for at in obj.atomList() if at.type.name.strip().lower() in ['sulphur','sulfur']]
2275 for s in sulphurs:
2276 neighbours = s.bondedTo()
2277 hydrogens = [neigh for neigh in neighbours if neigh.type.name.strip().lower() == 'hydrogen']
2278
2279 if len(hydrogens) == 1:
2280 selection.extend([o.index] + [h.index for h in hydrogens])
2281
2282 elif (len(hydrogens) == 2) and (obj.numberOfAtoms() == 3):
2283 selection.extend([s.index] + [h.index for h in hydrogens])
2284
2285 elif selectionKeyword == 'resname':
2286 for v in selectionValue:
2287 if v == '*':
2288 for obj in selectedObjects:
2289 selection.extend([at.index for at in obj.atomList()])
2290 else:
2291 for obj in selectedObjects:
2292 for res in obj.residues():
2293 if res.fullName().strip().lower() == v:
2294 selection.extend([at.index for at in res.atomList()])
2295
2296 elif selectionKeyword == 'restype':
2297 for v in selectionValue:
2298 if v == '*':
2299 for obj in selectedObjects:
2300 selection.extend([at.index for at in obj.atomList()])
2301 else:
2302 for obj in selectedObjects:
2303 for res in obj.residues():
2304 if res.symbol.strip().lower() == v:
2305 selection.extend([at.index for at in res.atomList()])
2306
2307 elif selectionKeyword == 'resclass':
2308 for v in selectionValue:
2309 for obj in selectedObjects:
2310 for res in obj.residues():
2311 if res.symbol.strip().lower() in residusChemFamily[v]:
2312 selection.extend([at.index for at in res.atomList()])
2313
2314 elif selectionKeyword == 'chainname':
2315 for v in selectionValue:
2316 if v == '*':
2317 for obj in selectedObjects:
2318 selection.extend([at.index for at in obj.atomList()])
2319 else:
2320 for obj in selectedObjects:
2321 for chain in obj:
2322 if chain.name.strip().lower() == v:
2323 selection.extend([at.index for at in chain.atomList()])
2324
2325 elif selectionKeyword == 'misc':
2326 for v in selectionValue:
2327 if v == 'sidechains':
2328 for obj in selectedObjects:
2329 selection.extend([at.index for at in obj.sidechains().atomList()])
2330
2331 elif v == 'backbone':
2332 for obj in selectedObjects:
2333 selection.extend([at.index for at in obj.backbone().atomList()])
2334
2335 return selection
2336
2337 -def setUniverseContents(universe):
2338 """Sets the contents of each object found in the universe.
2339
2340 @param universe: the MMTK universe to look in.
2341 @type universe: a instance of MMTK.Universe.
2342 """
2343
2344 if hasattr(universe, 'nmoldyncontents'):
2345 return
2346
2347 isotope = {}
2348 universe.indexToAtoms = {}
2349
2350 for obj in universe.objectList():
2351 if not isinstance(obj,(Atom, AtomCluster, Molecule, NucleotideChain, PeptideChain, Protein)):
2352 continue
2353
2354 for at in obj.atomList():
2355 try:
2356 universe.indexToAtoms[at.index] = at
2357 except:
2358 pass
2359
2360 try:
2361
2362 if obj.type.name:
2363 obj.nmoldynname = obj.type.name
2364 else:
2365 raise
2366
2367 except:
2368
2369 if isinstance(obj, Atom):
2370 obj.nmoldynname = 'A'
2371 elif isinstance(obj, AtomCluster):
2372 obj.nmoldynname = 'AC'
2373 elif isinstance(obj, Molecule):
2374 obj.nmoldynname = 'M'
2375 elif isinstance(obj, NucleotideChain):
2376 obj.nmoldynname = 'NC'
2377 elif isinstance(obj, PeptideChain):
2378 obj.nmoldynname = 'PC'
2379 elif isinstance(obj, Protein):
2380 obj.nmoldynname = 'P'
2381
2382 obj.nmoldynname += str(obj.numberOfAtoms())
2383
2384
2385 obj.nmoldynnamelist = sorted([at.name for at in obj.atomList()])
2386 if not isotope.has_key(obj.nmoldynname):
2387 isotope[obj.nmoldynname] = [obj.nmoldynnamelist]
2388
2389 else:
2390 if not obj.nmoldynnamelist in isotope[obj.nmoldynname]:
2391 isotope[obj.nmoldynname].append(obj.nmoldynnamelist)
2392
2393 universe.nmoldyncontents = {}
2394
2395
2396 for obj in universe.objectList():
2397
2398 if not isinstance(obj,(Atom, AtomCluster, Molecule, NucleotideChain, PeptideChain, Protein)):
2399 continue
2400
2401 if len(isotope[obj.nmoldynname]) > 1:
2402 LogMessage('info','Isotope found for %s' % obj.nmoldynname,['file','console'])
2403 iso = str(isotope[obj.nmoldynname].index(obj.nmoldynnamelist) + 1)
2404 obj.nmoldynname += '_iso' + iso
2405 obj.writeToFile(os.path.join(PREFERENCES.outputfile_path,obj.nmoldynname + '.pdb'), format = 'pdb')
2406
2407 delattr(obj, 'nmoldynnamelist')
2408
2409 if not obj.name:
2410 obj.name = obj.nmoldynname
2411
2412
2413
2414 if universe.nmoldyncontents.has_key(obj.nmoldynname):
2415 universe.nmoldyncontents[obj.nmoldynname]['name'].add(obj.name)
2416 universe.nmoldyncontents[obj.nmoldynname]['number'] += 1
2417
2418 else:
2419 universe.nmoldyncontents[obj.nmoldynname] = {'name' : set(['*', obj.name]) , 'number' : 1}
2420
2421
2422
2423 if isinstance(obj, Atom):
2424 universe.nmoldyncontents[obj.nmoldynname]['objectclass'] = 'Atom'
2425 universe.nmoldyncontents[obj.nmoldynname]['groupinglevel'] = [0, ['atom',]]
2426
2427
2428
2429
2430
2431 elif isinstance(obj, AtomCluster):
2432 universe.nmoldyncontents[obj.nmoldynname]['objectclass'] = 'AtomCluster'
2433 universe.nmoldyncontents[obj.nmoldynname]['groupinglevel'] = [1, ['atom','cluster']]
2434 universe.nmoldyncontents[obj.nmoldynname]['atomname'] = set(['*'])
2435 universe.nmoldyncontents[obj.nmoldynname]['atomelement'] = set(['*'])
2436 for atom in obj.atomList():
2437 universe.nmoldyncontents[obj.nmoldynname]['atomname'].add(atom.fullName())
2438 universe.nmoldyncontents[obj.nmoldynname]['atomelement'].add(atom.type.name)
2439
2440 universe.nmoldyncontents[obj.nmoldynname]['atomname'] = sorted(universe.nmoldyncontents[obj.nmoldynname]['atomname'])
2441 universe.nmoldyncontents[obj.nmoldynname]['atomelement'] = sorted(universe.nmoldyncontents[obj.nmoldynname]['atomelement'])
2442
2443
2444
2445
2446
2447
2448 elif isinstance(obj, Molecule):
2449 universe.nmoldyncontents[obj.nmoldynname]['objectclass'] = 'Molecule'
2450 universe.nmoldyncontents[obj.nmoldynname]['groupinglevel'] = [5, ['atom','amine', 'hydroxy', 'methyl', 'thiol', 'molecule']]
2451 universe.nmoldyncontents[obj.nmoldynname]['atomname'] = set(['*'])
2452 universe.nmoldyncontents[obj.nmoldynname]['atomelement'] = set(['*'])
2453 universe.nmoldyncontents[obj.nmoldynname]['chemfragment'] = ['amine', 'hydroxy', 'methyl', 'thiol']
2454
2455 for atom in obj.atomList():
2456 universe.nmoldyncontents[obj.nmoldynname]['atomname'].add(atom.fullName())
2457 universe.nmoldyncontents[obj.nmoldynname]['atomelement'].add(atom.type.name)
2458
2459 universe.nmoldyncontents[obj.nmoldynname]['atomname'] = sorted(universe.nmoldyncontents[obj.nmoldynname]['atomname'])
2460 universe.nmoldyncontents[obj.nmoldynname]['atomelement'] = sorted(universe.nmoldyncontents[obj.nmoldynname]['atomelement'])
2461
2462
2463
2464
2465
2466
2467
2468
2469
2470 elif isinstance(obj, NucleotideChain):
2471 universe.nmoldyncontents[obj.nmoldynname]['objectclass'] = 'NucleotideChain'
2472 universe.nmoldyncontents[obj.nmoldynname]['groupinglevel'] = [2, ['atom', 'amine', 'residue', 'nucleicacid']]
2473 universe.nmoldyncontents[obj.nmoldynname]['atomname'] = ['*']
2474 universe.nmoldyncontents[obj.nmoldynname]['atomtype'] = set(['*'])
2475 universe.nmoldyncontents[obj.nmoldynname]['atomelement'] = set(['*'])
2476 universe.nmoldyncontents[obj.nmoldynname]['nuclname'] = ['*']
2477 universe.nmoldyncontents[obj.nmoldynname]['nucltype'] = set(['*'])
2478 universe.nmoldyncontents[obj.nmoldynname]['misc'] = ['backbone', 'bases']
2479
2480
2481 for nucl in obj.residues():
2482 universe.nmoldyncontents[obj.nmoldynname]['nuclname'].append(nucl.fullName())
2483 universe.nmoldyncontents[obj.nmoldynname]['nucltype'].add(nucl.symbol)
2484
2485
2486 for atom in nucl.atomList():
2487 universe.nmoldyncontents[obj.nmoldynname]['atomname'].append(atom.fullName())
2488 universe.nmoldyncontents[obj.nmoldynname]['atomtype'].add(atom.name)
2489 universe.nmoldyncontents[obj.nmoldynname]['atomelement'].add(atom.type.name)
2490
2491 universe.nmoldyncontents[obj.nmoldynname]['atomtype'] = sorted(universe.nmoldyncontents[obj.nmoldynname]['atomtype'])
2492 universe.nmoldyncontents[obj.nmoldynname]['atomelement'] = sorted(universe.nmoldyncontents[obj.nmoldynname]['atomelement'])
2493
2494 universe.nmoldyncontents[obj.nmoldynname]['nucltype'] = sorted(universe.nmoldyncontents[obj.nmoldynname]['nucltype'])
2495
2496
2497
2498
2499
2500
2501
2502
2503
2504
2505
2506 elif isinstance(obj, PeptideChain):
2507 universe.nmoldyncontents[obj.nmoldynname]['objectclass'] = 'PeptideChain'
2508 universe.nmoldyncontents[obj.nmoldynname]['groupinglevel'] = [5, ['atom', 'amine', 'hydroxy', 'methyl', 'thiol', 'residue','chain']]
2509 universe.nmoldyncontents[obj.nmoldynname]['atomname'] = ['*']
2510 universe.nmoldyncontents[obj.nmoldynname]['atomtype'] = set(['*'])
2511 universe.nmoldyncontents[obj.nmoldynname]['atomelement'] = set(['*'])
2512 universe.nmoldyncontents[obj.nmoldynname]['chemfragment'] = ['amine', 'c_alphas', 'hydroxy', 'methyl', 'thiol']
2513 universe.nmoldyncontents[obj.nmoldynname]['resname'] = ['*']
2514 universe.nmoldyncontents[obj.nmoldynname]['restype'] = set(['*'])
2515 universe.nmoldyncontents[obj.nmoldynname]['resclass'] = sorted(residusChemFamily.keys())
2516 universe.nmoldyncontents[obj.nmoldynname]['misc'] = ['backbone', 'sidechains']
2517
2518
2519 for residue in obj.residues():
2520 universe.nmoldyncontents[obj.nmoldynname]['resname'].append(residue.fullName())
2521 universe.nmoldyncontents[obj.nmoldynname]['restype'].add(residue.symbol)
2522
2523
2524 for atom in residue.atomList():
2525 universe.nmoldyncontents[obj.nmoldynname]['atomname'].append(atom.fullName())
2526 universe.nmoldyncontents[obj.nmoldynname]['atomtype'].add(atom.name)
2527 universe.nmoldyncontents[obj.nmoldynname]['atomelement'].add(atom.type.name)
2528
2529 universe.nmoldyncontents[obj.nmoldynname]['atomtype'] = sorted(universe.nmoldyncontents[obj.nmoldynname]['atomtype'])
2530 universe.nmoldyncontents[obj.nmoldynname]['atomelement'] = sorted(universe.nmoldyncontents[obj.nmoldynname]['atomelement'])
2531
2532 universe.nmoldyncontents[obj.nmoldynname]['restype'] = sorted(universe.nmoldyncontents[obj.nmoldynname]['restype'])
2533
2534
2535
2536
2537
2538
2539
2540
2541
2542
2543
2544
2545 elif isinstance(obj, Protein):
2546 universe.nmoldyncontents[obj.nmoldynname]['objectclass'] = 'Protein'
2547 universe.nmoldyncontents[obj.nmoldynname]['groupinglevel'] = [5, ['atom', 'amine', 'hydroxy', 'methyl', 'thiol', 'residue','chain', 'protein']]
2548 universe.nmoldyncontents[obj.nmoldynname]['atomname'] = ['*']
2549 universe.nmoldyncontents[obj.nmoldynname]['atomtype'] = set(['*'])
2550 universe.nmoldyncontents[obj.nmoldynname]['atomelement'] = set(['*'])
2551 universe.nmoldyncontents[obj.nmoldynname]['chemfragment'] = ['amine', 'c_alphas', 'hydroxy', 'methyl', 'thiol']
2552 universe.nmoldyncontents[obj.nmoldynname]['resname'] = ['*']
2553 universe.nmoldyncontents[obj.nmoldynname]['restype'] = set(['*'])
2554 universe.nmoldyncontents[obj.nmoldynname]['resclass'] = sorted(residusChemFamily.keys())
2555 universe.nmoldyncontents[obj.nmoldynname]['chainname'] = ['*']
2556 universe.nmoldyncontents[obj.nmoldynname]['misc'] = ['backbone', 'sidechains']
2557
2558
2559 for chain in obj:
2560
2561 universe.nmoldyncontents[obj.nmoldynname]['chainname'].append(chain.fullName())
2562
2563
2564 for residue in chain.residues():
2565 universe.nmoldyncontents[obj.nmoldynname]['resname'].append(residue.fullName())
2566 universe.nmoldyncontents[obj.nmoldynname]['restype'].add(residue.symbol)
2567
2568
2569 for atom in residue.atomList():
2570 universe.nmoldyncontents[obj.nmoldynname]['atomname'].append(atom.fullName())
2571 universe.nmoldyncontents[obj.nmoldynname]['atomtype'].add(atom.name)
2572 universe.nmoldyncontents[obj.nmoldynname]['atomelement'].add(atom.type.name)
2573
2574 universe.nmoldyncontents[obj.nmoldynname]['atomtype'] = sorted(universe.nmoldyncontents[obj.nmoldynname]['atomtype'])
2575 universe.nmoldyncontents[obj.nmoldynname]['atomelement'] = sorted(universe.nmoldyncontents[obj.nmoldynname]['atomelement'])
2576
2577 universe.nmoldyncontents[obj.nmoldynname]['restype'] = sorted(universe.nmoldyncontents[obj.nmoldynname]['restype'])
2578
2579
2580 for objTypeName in universe.nmoldyncontents.keys():
2581 universe.nmoldyncontents[objTypeName]['name'] = sorted(universe.nmoldyncontents[objTypeName]['name'])
2582
2583
2584 universe.nmoldyncontents['*'] = {'atomelement' : set(['*']),\
2585 'number' : universe.numberOfAtoms(),\
2586 'objectclass' : 'AllClass',\
2587 'groupinglevel' : [0,['default',]]}
2588 for atom in universe.atomList():
2589 universe.nmoldyncontents['*']['atomelement'].add(atom.type.name)
2590 universe.nmoldyncontents['*']['atomelement'] = sorted(universe.nmoldyncontents['*']['atomelement'])
2591
2593 """Generates a set of QVectors within a given shell.
2594 """
2595 - def __init__(self, universe, generator, qRadii, dq, qVectorsPerShell, qVectorsDirection = None):
2596 """The constructor.
2597
2598 @param universe: the MMTK universe used to define the reciprocal space.
2599 @type universe: a MMTK.Universe subclass object
2600
2601 @param generator: a string being one of '3d isotropic', '2d isotropic' or 'anistropic' the way the q-vectors
2602 should be generated.
2603 @type generator: string
2604
2605 @param qRadii: a list of floats specifying the radii of the shell in which the q vectors have to be
2606 generated.
2607 @type qRadii: list
2608
2609 @param dq: a float specifying the width of a qhsell defined as [|qRadius| - dq/2,|qRadius|+dq/2].
2610 @type dq: float
2611
2612 @param qVectorsPerShell: an integer specifying the number of q-vectors to generate for each shell.
2613 @type qVectorsPerShell: integer
2614
2615 @param qVectorsDirection: a list of Scientific.Geometry.Vector objects specifying the directions along which
2616 the q-vectors should be generated. If None, the q-vectors generation will be isotropic.
2617 @type qVectorsDirection: list
2618 """
2619
2620 self.universe = universe
2621 self.generator = generator
2622 self.dq = dq
2623 self.qVectorsPerShell = qVectorsPerShell
2624 self.qVectorsDirection = qVectorsDirection
2625
2626 if self.generator == '3d isotropic':
2627 self.qVectorsDirection = None
2628
2629 else:
2630
2631 if self.qVectorsDirection is None:
2632 raise Error('Directions must be given for %s q vectors generation.' % self.generator)
2633
2634 if self.generator == '2d isotropic':
2635 if len(self.qVectorsDirection) != 2:
2636 raise Error('Two vectors must be given to generate q vectors in a plane.')
2637
2638 v1, v2 = self.qVectorsDirection
2639 if abs(v1.angle(v2)) < 1.e-10:
2640 raise Error('The two vectors are colinear.')
2641
2642 self.qVectorsDirection = qVectorsDirection
2643
2644 self.qRadii = []
2645 self.qMin = []
2646 self.qMax = []
2647
2648 for qRadius in sorted(qRadii):
2649
2650 if qRadius <= 0.0:
2651 LogMessage('warning','No Q vector could be generated for shell with radius Q = %s.' % qRadius,['file','console'])
2652 continue
2653 else:
2654 self.qRadii.append(qRadius)
2655 self.qMin.append(max(qRadius - 0.5*self.dq,0.0))
2656 self.qMax.append(qRadius + 0.5*self.dq)
2657
2658
2659 self.reciprocalBasis = universe.reciprocalBasisVectors()
2660
2661
2662 if self.reciprocalBasis is None:
2663 self.directMatrix = None
2664 self.reciprocalCellVolume = None
2665
2666
2667 else:
2668 self.reciprocalBasis = [2.0*N.pi*v for v in self.reciprocalBasis]
2669 self.directMatrix = N.array([N.array(v) for v in self.universe.basisVectors()])
2670 self.directMatrix = Tensor(self.directMatrix)/(2.*N.pi)
2671 self.reciprocalCellVolume = self.reciprocalBasis[0]*(self.reciprocalBasis[1].cross(self.reciprocalBasis[2]))
2672
2673 self.qVectors = {}
2674
2675
2676 for index in range(len(self.qRadii)):
2677 qRadius = self.qRadii[index]
2678 self.qVectors[qRadius] = []
2679 qMin = self.qMin[index]
2680 qMax = self.qMax[index]
2681
2682 if self.generator in ['2d isotropic', '3d isotropic']:
2683 self.__isotropicGeneration(qRadius, qMin, qMax)
2684
2685 elif self.generator == '2d isotropic':
2686 self.__planarIsotropicGeneration(qRadius, qMin, qMax)
2687
2688 elif self.generator == 'anisotropic':
2689 self.__anisotropicGeneration(qRadius, qMin, qMax)
2690
2691 if len(self.qVectors[qRadius]) == 0:
2692 LogMessage('warning','No Q vector could be generated for shell with radius Q = %s.' % qRadius,['file','console'])
2693 del self.qVectors[qRadius]
2694 continue
2695
2696 if not self.qVectors:
2697 raise Error('No Q Vectors generated. Please check your settings.')
2698
2699 self.statistics = N.zeros((len(self.qVectors), 8), typecode = N.Int)
2700 comp = 0
2701 for qRadius in sorted(self.qVectors.keys()):
2702 qVectors = self.qVectors[qRadius]
2703 for qVect in qVectors:
2704 if qVect[0] > 0:
2705 if qVect[1] > 0:
2706
2707 if qVect[2] > 0:
2708 self.statistics[comp, 0] += 1
2709
2710 else:
2711 self.statistics[comp, 1] += 1
2712 else:
2713
2714 if qVect[2] > 0:
2715 self.statistics[comp, 2] += 1
2716
2717 else:
2718 self.statistics[comp, 3] += 1
2719 else:
2720 if qVect[1] > 0:
2721
2722 if qVect[2] > 0:
2723 self.statistics[comp, 4] += 1
2724
2725 else:
2726 self.statistics[comp, 5] += 1
2727 else:
2728
2729 if qVect[2] > 0:
2730 self.statistics[comp, 6] += 1
2731
2732 else:
2733 self.statistics[comp, 7] += 1
2734 comp += 1
2735
2737 if self.reciprocalCellVolume is None:
2738 self.__makeRandomQList()
2739
2740 else:
2741 reciprocalShellVolume = 4.0*N.pi*qRadius**2 * (qMax - qMin)
2742 nCells = reciprocalShellVolume/self.reciprocalCellVolume
2743 if nCells > 500:
2744 LogMessage('info', 'Random Q vectors generation for shell with radius Q = %s.' % qRadius, ['file','console'])
2745 self.__makeRandomQList(qRadius, qMin, qMax)
2746 else:
2747 LogMessage('info', 'Explicit Q vectors generation for shell with radius Q = %s.' % qRadius, ['file','console'])
2748 self.__makeExplicitQList(qRadius, qMin, qMax)
2749
2751 LogMessage('info', 'Explicit Q vectors generation for shell with radius Q = %s.' % qRadius, ['file','console'])
2752 self.__makeExplicitQList(qRadius, qMin, qMax)
2753
2755
2756
2757
2758 redundancyCounter = 0
2759 while len(self.qVectors[qRadius]) < self.qVectorsPerShell:
2760
2761
2762 q = randomVector(self.qVectorsDirection)*uniform(qMin, qMax)
2763 if self.reciprocalBasis is not None:
2764
2765
2766 q = [round(v) for v in N.array(self.directMatrix*q)]
2767
2768 q = q[0]*self.reciprocalBasis[0] + q[1]*self.reciprocalBasis[1] + q[2]*self.reciprocalBasis[2]
2769
2770
2771
2772 if qMin < q.length() <= qMax:
2773
2774 if q in self.qVectors[qRadius]:
2775 redundancyCounter += 1
2776
2777 if redundancyCounter == 5000:
2778 return
2779 continue
2780 self.qVectors[qRadius].append(q)
2781
2782 redundancyCounter = 0
2783 else:
2784 redundancyCounter += 1
2785
2786 if redundancyCounter == 5000:
2787 return
2788
2790
2791
2792 if self.generator == '3d isotropic':
2793 upperLimit = [int(qMax/v.length()) for v in self.reciprocalBasis]
2794
2795 for l in range(-upperLimit[0], upperLimit[0]+1):
2796 for m in range(-upperLimit[1], upperLimit[1]+1):
2797 for n in range(-upperLimit[2], upperLimit[2]+1):
2798 if l == m == n == 0:
2799 continue
2800
2801 q = l*self.reciprocalBasis[0] + m*self.reciprocalBasis[1] + n*self.reciprocalBasis[2]
2802
2803 if qMin < q.length() <= qMax:
2804 if q in self.qVectors[qRadius]:
2805 continue
2806 self.qVectors[qRadius].append(q)
2807
2808 elif self.generator == '2d isotropic':
2809 upperLimit = [int(qMax/v.length()) for v in self.qVectorsDirection]
2810 for l in range(-upperLimit[0], upperLimit[0]+1):
2811 for m in range(-upperLimit[1], upperLimit[1]+1):
2812 if l == m == 0:
2813 continue
2814
2815
2816 q = l*self.qVectorsDirection[0] + m*self.qVectorsDirection[1]
2817
2818
2819
2820 q = [round(v) for v in N.array(self.directMatrix*q)]
2821
2822 q = q[0]*self.reciprocalBasis[0] + q[1]*self.reciprocalBasis[1] + q[2]*self.reciprocalBasis[2]
2823
2824 if qMin < q.length() <= qMax:
2825 if q in self.qVectors[qRadius]:
2826 continue
2827 self.qVectors[qRadius].append(q)
2828
2829 elif self.generator == 'anisotropic':
2830
2831 for qVect in self.qVectorsDirection:
2832 upperLimit = int(qMax/qVect.length())
2833 for l in range(-upperLimit, upperLimit + 1):
2834 if l == 0:
2835 continue
2836
2837
2838
2839 q = l*qVect
2840
2841
2842
2843 q = [round(v) for v in N.array(self.directMatrix*q)]
2844
2845 q = q[0]*self.reciprocalBasis[0] + q[1]*self.reciprocalBasis[1] + q[2]*self.reciprocalBasis[2]
2846
2847 if qMin < q.length() <= qMax:
2848 if q in self.qVectors[qRadius]:
2849 continue
2850 self.qVectors[qRadius].append(q)
2851
2852
2853
2854 if len(self.qVectors[qRadius]) > self.qVectorsPerShell:
2855 LogMessage('warning', '%s qVectors generated for shell with radius Q = %s. Will be reduced to %s.' % (len(self.qVectors[qRadius]), qRadius, self.qVectorsPerShell), ['file','console'])
2856 self.qVectors[qRadius] = sample(self.qVectors[qRadius], self.qVectorsPerShell)
2857