diff --git a/main.py b/main.py index 2e2e9c9..01dfd27 100755 --- a/main.py +++ b/main.py @@ -233,20 +233,27 @@ class MayaviPhononWindow(QtGui.QMainWindow): def __init__(self, crystal_structure,data_dictionary, parent=None): super(MayaviPhononWindow, self).__init__(parent) + self.resize(1000, 600) + self.data_dictionary = data_dictionary self.setWindowTitle('Phonon visualization') self.main_widget = QtGui.QWidget(self) self.setCentralWidget(self.main_widget) - layout = QtGui.QHBoxLayout(self.main_widget) + main_layout = QtGui.QVBoxLayout(self.main_widget) + + self.sub_widget = QtGui.QWidget(self.main_widget) + main_layout.addWidget(self.sub_widget) + + layout = QtGui.QHBoxLayout(self.sub_widget) layout.setContentsMargins(0, 0, 0, 0) layout.setSpacing(0) self.visualization = PhononVisualization(crystal_structure) - self.ui = self.visualization.edit_traits(parent=self.main_widget, + self.ui = self.visualization.edit_traits(parent=self.sub_widget, kind='subpanel').control layout.addWidget(self.ui) - self.ui.setParent(self.main_widget) + self.ui.setParent(self.sub_widget) self.treeview = QtGui.QTreeWidget(parent=self) self.treeview.setMaximumWidth(200) @@ -266,6 +273,21 @@ def __init__(self, crystal_structure,data_dictionary, parent=None): self.k_treeview.itemSelectionChanged.connect(self.handle_item_changed) layout.addWidget(self.k_treeview) + self.info_widget = QtGui.QWidget(self.main_widget) + main_layout.addWidget(self.info_widget) + + info_layout = QtGui.QHBoxLayout(self.info_widget) + info_layout.setAlignment(QtCore.Qt.AlignLeft) + + infos = ['frequency'] + self.info_labels = {} + + for info in infos: + new_label = LabeledLabel(info.title()+':') + info_layout.addWidget(new_label) + self.info_labels[info] = new_label + + def update_plot(self, keep_view=False): self.visualization.update_plot(keep_view=keep_view) @@ -312,6 +334,9 @@ def handle_item_changed(self): mode,k = self.get_mode_and_k() self.visualization.plot_phonons(mode,k) + self.update_infos(phonon_eigenvectors,mode,k) + + def add_result_key(self, title): item = QtGui.QTreeWidgetItem(self.treeview.invisibleRootItem(), [title]) return item @@ -370,6 +395,16 @@ def get_mode_and_k(self): return mode,k + def update_infos(self,eigs,mode,k): + freqs = eigs.frequencies + freq = freqs[mode,k] + data = {'frequency':'{0:1.1f} cm^-1'.format(freq)} + + for label,widget in self.info_labels.items(): + if label in data.keys(): + widget.text(data[label]) + + class EntryWithLabel(QtGui.QWidget): def __init__(self, parent, label, value=None, width_text=200, width_label=90): @@ -2747,7 +2782,7 @@ def __init__(self, parent=None, *args, **kwargs): if DEBUG: if sys.platform in ['linux', 'linux2']: - project_directory = r"/home/jannick/OpenDFT_projects/diamond3/" + project_directory = r"/home/jannick/OpenDFT_projects/CH4_exciting/" # project_directory = r"/home/jannick/exciting_cluster/GaN" else: project_directory = r'D:\OpenDFT_projects\test' diff --git a/nwchem_handler.py b/nwchem_handler.py index 8e08530..3c6e424 100644 --- a/nwchem_handler.py +++ b/nwchem_handler.py @@ -63,7 +63,7 @@ def __init__(self): self._timestamp_tasks = {} self.supported_methods = sst.ComputationalMethods( - ['non-periodic', 'scf', 'relax']) + ['non-periodic', 'scf', 'relax','phonons']) self.project_directory = None self.input_filename = 'scf.in' @@ -74,10 +74,11 @@ def __init__(self): self.custom_command = '' self.custom_command_active = False self.dft_installation_folder = self.find_engine_folder() - self.scf_options = convert_to_ordered( - {'basis': 'cc-pvdz', 'method': 'scf', 'spin state': 'singlet', 'spin restriction': 'RHF'}) + self.scf_options = convert_to_ordered({'basis': 'cc-pvdz', 'method': 'scf', + 'spin state': 'singlet', 'spin restriction': 'RHF', 'xc': 'pbe0','charge':'0'}) self.scf_options_tooltip = {} + self.general_options = {'title': 'title'} self.bs_options = {} self.relax_options = {} @@ -157,7 +158,7 @@ def start_gw(self, crystal_structure, band_structure_points=None, blocking=False """ raise NotImplementedError - def start_phonon(self, crystal_structure, band_structure_points): + def start_phonon(self, crystal_structure, band_structure_points,blocking=False): """This method starts a phonon bandstructure calculation in a subprocess. The configuration is stored in phonons_options. Args: @@ -171,9 +172,12 @@ def start_phonon(self, crystal_structure, band_structure_points): Returns: - None """ - raise NotImplementedError + file = self._make_input_file() + self._add_scf_to_file(file, crystal_structure,calculation='freq') + file.close() + self._start_engine(blocking=blocking) - def start_relax(self, crystal_structure): + def start_relax(self, crystal_structure,blocking=False): """This method starts a structure relaxation calculation in a subprocess. The configuration is stored in relax_options. Args: @@ -185,7 +189,7 @@ def start_relax(self, crystal_structure): file = self._make_input_file() self._add_scf_to_file(file, crystal_structure, calculation='optimize') file.close() - self._start_engine() + self._start_engine(blocking=blocking) def load_relax_structure(self): """This method loads the result of a relaxation calculation, which is a molecular or crystal structure. @@ -239,43 +243,45 @@ def read_scf_status(self): return None info_text = f.read() f.close() + try: + scf_energy_list = [] - scf_energy_list = [] - - # this is for future use - matches = re.findall('SCF calculation type: \w*', info_text) - if len(matches) == 0: - calculation_type = 'scf' - task = 'scf' - else: - calculation_type = matches[0].split(':')[1].strip().lower() - task = '' - - hits = [i for i, l in enumerate(info_text.split('\n')) if - l.strip().startswith('iter') or l.strip().startswith('convergence')] - if len(hits) > 0: - hit = hits[0] - header_line = info_text.split('\n')[hit].split() - energy_index = header_line.index('energy') - - for line in info_text.split('\n')[hit + 2:]: - sline = line.split(' ') - sline = [x for x in sline if len(x) > 0] - if len(sline) <= energy_index: - break - scf_energy_list.append(float(sline[energy_index].split()[0])) - - # find only final result - matches = re.findall(r"Total \w\w\w energy[\s\t]*=[\s\t]*[-+]?\d*\.\d+", info_text) - for match in matches: - ms = match.split('=') - scf_energy_list.append(float(ms[1])) - - res = np.array(list(zip(range(1, len(scf_energy_list) + 1), scf_energy_list))) - - if len(res) == 0: + # this is for future use + matches = re.findall('SCF calculation type: \w*', info_text) + if len(matches) == 0: + calculation_type = 'scf' + task = 'scf' + else: + calculation_type = matches[0].split(':')[1].strip().lower() + task = '' + + hits = [i for i, l in enumerate(info_text.split('\n')) if + l.strip().startswith('iter') or l.strip().startswith('convergence')] + if len(hits) > 0: + hit = hits[0] + header_line = info_text.split('\n')[hit].split() + energy_index = header_line.index('energy') + + for line in info_text.split('\n')[hit + 2:]: + sline = line.split(' ') + sline = [x for x in sline if len(x) > 0] + if len(sline) <= energy_index: + break + scf_energy_list.append(float(sline[energy_index].split()[0])) + + # find only final result + matches = re.findall(r"Total \w\w\w energy[\s\t]*=[\s\t]*[-+]?\d*\.\d+", info_text) + for match in matches: + ms = match.split('=') + scf_energy_list.append(float(ms[1])) + + res = np.array(list(zip(range(1, len(scf_energy_list) + 1), scf_energy_list))) + + if len(res) == 0: + return None + return res + except Exception as e: return None - return res def read_bandstructure(self, special_k_points=None): """This method reads the result of a electronic band structure calculation. @@ -381,13 +387,60 @@ def read_gw_bandstructure(self, filename='BAND-QP.OUT'): """ raise NotImplementedError - def read_phonon_bandstructure(self): + def read_phonon_bandstructure(self,special_k_points=None, structure=None): """This method reads the result of a phonon band structure calculation. Returns: - band_structure: A BandStructure object with the latest phonon band structure result found. """ - raise NotImplementedError + try: + f = open(self.project_directory + self.working_dirctory + self.info_file, 'r') + except IOError: + return None + text = f.readlines() + + frequencies = [] + eigenvectors = [] + + freq_found = False + eig_lines = [] + + eigenvector_array = None + + for line in text: + if line.strip().startswith('P.Frequency'): + freq_strings = line.split()[1:] + freqs = [float(x) for x in freq_strings] + frequencies.extend(freqs) + freq_found = True + continue + + + if freq_found: + if len(line.strip()) == 0: + continue + else: + eig_line = line.strip().split() + if len(eig_line) == len(freqs)+1: + eig_line_float = [float(x) for x in eig_line[1:]] + eig_lines.append(eig_line_float) + else: + eig_array = np.array(eig_lines).T + if eigenvector_array is None: + eigenvector_array = eig_array + else: + eigenvector_array = np.append(eigenvector_array,eig_array,axis=0) + freqs = [] + eig_lines = [] + freq_found = False + k_points = np.array([0,0,0]) + k_points = k_points[np.newaxis,:] + freqs_conv = np.array(frequencies)[:,np.newaxis] + eigs_conv = eigenvector_array[:,np.newaxis,:] + phonon_eigenvectors = sst.PhononEigenvectors(freqs_conv,eigs_conv,k_points) + vib = sst.VibrationalStructure(np.array(frequencies)) + vib.phonon_eigenvectors = phonon_eigenvectors + return vib def read_optical_spectrum(self): """This method reads the result of a optical spectrum calculation. @@ -525,9 +578,12 @@ def _add_scf_to_file(self, file, crystal_structure, calculation='scf', band_poin file.write('title ' + '"' + self.general_options['title'] + '"\n') self._add_geometry(file, crystal_structure) self._add_basis(file, crystal_structure) - self._add_scf_field_to_file(file) + if self.scf_options['method'] =='scf': + self._add_scf_field_to_file(file) + elif self.scf_options['method'] == 'dft': + self._add_dft_field_to_file(file,crystal_structure) file.write('task ' + self.scf_options['method']) - if calculation == 'optimize': + if calculation in ['optimize','freq']: file.write(' ' + calculation) file.write('\n') @@ -566,6 +622,7 @@ def _add_geometry(self, file, crystal_structure, auto=False): file.write(' ' + p_table[atom[3]].lower()) file.write(" {0:1.9f} {1:1.9f} {2:1.9f}\n".format(*atom[:3])) file.write('end\n') + file.write('charge '+self.scf_options['charge']+'\n') def _add_basis(self, file, crystal_structure): file.write('basis\n') @@ -583,6 +640,17 @@ def _add_scf_field_to_file(self, file, input=False): file.write('end\n') + def _add_dft_field_to_file(self, file,crystal_structure): + n_electrons = int(crystal_structure.atoms[:,3].sum()) + mult_dic = {'singlet':1, 'doublet':2, 'triplet':3,'quartet':4,'quintet':5,'sextet':6,'septet':7,'octet':8} + + file.write('dft\n') + file.write(' xc ' + self.scf_options['xc'] + '\n') + if n_electrons%2 == 1: + file.write(' ODFT' + '\n') + file.write(' MULT {}'.format(mult_dic[self.scf_options['spin state']])+'\n') + file.write('end\n') + def _add_dplot_to_file(self, file, crystal_structure, limits=None, orbital=None, grid=(100, 100, 100)): if limits is None: @@ -618,11 +686,14 @@ def _add_dplot_to_file(self, file, crystal_structure, limits=None, orbital=None, crystal_structure = sst.MolecularStructure(atoms) handler = Handler() - handler.project_directory = "/home/jannick/OpenDFT_projects/nwchem3" + handler.project_directory = "/home/jannick/OpenDFT_projects/CH4_exciting" + + eigs = handler.read_phonon_bandstructure() + # handler.start_ground_state(crystal_structure) # handler.read_scf_status() - energies = handler.read_scf_status() - print(energies) + # energies = handler.read_scf_status() + # print(energies) # handler.calculate_electron_density(crystal_structure) # while handler.is_engine_running(): # time.sleep(0.1) diff --git a/solid_state_tools.py b/solid_state_tools.py index 82e7900..ae5fccd 100755 --- a/solid_state_tools.py +++ b/solid_state_tools.py @@ -264,6 +264,12 @@ def _find_bandgap(self, bands): return bandgap, k_bandgap +class VibrationalStructure(object): + def __init__(self,frequencies): + self.frequencies = frequencies + self.engine_information = None + + class EnergyDiagram(object): def __init__(self, energies, labels, occupations=None): self.energies = energies diff --git a/visualization.py b/visualization.py index 099df2d..8888e46 100755 --- a/visualization.py +++ b/visualization.py @@ -7,7 +7,7 @@ os.environ['ETS_TOOLKIT'] = 'qt4' from pyface.qt import QtGui, QtCore -from traits.api import HasTraits, Instance, on_trait_change, Range, Bool, Button, Array +from traits.api import HasTraits, Instance, on_trait_change, Range, Bool, Button, Array, Float, Enum from traitsui.api import View, Item, Group, HGroup from mayavi import mlab @@ -1010,6 +1010,8 @@ def plot(self, band_structure, *args, **kwargs): self.plot_energy_diagram(band_structure) elif type(band_structure) is sst.BandStructure: self.plot_bandstructure(band_structure) + elif type(band_structure) is sst.VibrationalStructure: + self.plot_vibrational(band_structure) try: Emin = float(self.Emin_entry.get_text()) @@ -1121,6 +1123,24 @@ def plot_bandstructure(self, band_structure): self.ax.set_ylim(bottom=0) self.ax.set_title('Phonon bandstructure', fontsize=25, color=title_color) + def plot_vibrational(self,band_structure): + self.ax.format_coord = lambda x, y: r'ν = {0:1.1f} cm^-1'.format(y) + + self.ax.cla() + if self.dark_mode: + set_dark_mode_matplotlib(self.figure, self.ax, self.bg_color) + title_color = 'white' + else: + title_color = 'black' + + for i,freq in enumerate(band_structure.frequencies): + self.ax.plot([-.3, .3],[freq,freq]) + + self.ax.set_ylabel(r'Frequncy $\left[\mathrm{cm}^{-1}\right]$') + self.ax.set_xticks([]) + self.ax.set_xlim(-1, 1) + self.ax.set_title('Vibrational modes', fontsize=25, color=title_color) + def make_interactive_text(self, k_in, E, band_structure): bands = band_structure.bands k = band_structure.bands[0][:, 0] @@ -1349,9 +1369,12 @@ def plot(self, scf_data): class PhononVisualization(StructureVisualization): + arrow = Float(10) + colormap = Enum('viridis',colormap_list) + view = View(Item('scene', editor=SceneEditor(scene_class=MayaviScene), height=450, width=500, show_label=False), - Group('_','show_unitcell', 'show_bonds', 'show_atoms', orientation='horizontal'), + Group('_','show_unitcell', 'show_bonds', 'show_atoms','arrow','colormap', orientation='horizontal'), resizable=True, # We need this to resize with the parent widget ) @@ -1361,12 +1384,19 @@ def __init__(self,crystal_structure): super(PhononVisualization,self).__init__(crystal_structure) self.phonon_eigenvectors = None self.mayavi_phonons = None + self.last_plot = None + @on_trait_change('arrow,colormap') + def _arrow_changed_event(self): + if self.last_plot: + self.plot_phonons(*self.last_plot) def plot_phonons(self,n_mode,n_k): self.remove_phonons() if self.phonon_eigenvectors is None: return + self.last_plot = [n_mode,n_k] + repeat = [self.n_x, self.n_y, self.n_z] abs_coord_atoms = self.crystal_structure.calc_absolute_coordinates(repeat=repeat, edges=self.edge_atoms) @@ -1376,12 +1406,14 @@ def plot_phonons(self,n_mode,n_k): mode = mode_conv.reshape(self.crystal_structure.atoms.shape[0],3) # mode_unit = mode/np.abs(mode).max() self.mayavi_phonons = self.scene.mlab.quiver3d(abs_coord_atoms[:,0],abs_coord_atoms[:,1],abs_coord_atoms[:,2],mode[:,0],mode[:,1],mode[:,2],figure=self.scene.mayavi_scene, - scale_factor=10,line_width=3,reset_zoom=False,colormap='viridis') + scale_factor=self.arrow,line_width=3,reset_zoom=False,colormap=self.colormap) + def remove_phonons(self): if self.mayavi_phonons: self.mayavi_phonons.remove() self.mayavi_phonons = None + def set_dark_mode_matplotlib(f, ax, color): ax.tick_params(axis='x', colors='white') ax.tick_params(axis='y', colors='white')