From 37c9df3022a3d63be0b01b39be9f944af0fb36f5 Mon Sep 17 00:00:00 2001 From: abel Date: Fri, 13 Aug 2021 17:30:53 +0200 Subject: [PATCH] Added function for getting custom basis from BSE --- README.md | 4 +- examples/custom_basis.py | 17 +++++++ pyqchem/__init__.py | 2 +- pyqchem/basis.py | 98 ++++++++++++++++++++++++++++++++++++++-- 4 files changed, 115 insertions(+), 6 deletions(-) diff --git a/README.md b/README.md index b59a36b..bda8600 100644 --- a/README.md +++ b/README.md @@ -104,7 +104,7 @@ for i, mode in enumerate(parsed_data['modes']): ```python from pyqchem import QchemInput, Structure -from pyqchem.basis import get_basis_from_ccRepo +from pyqchem.basis import get_basis_from_BSE molecule = Structure(coordinates=[[0.0, 0.0, 0.0000], @@ -113,7 +113,7 @@ molecule = Structure(coordinates=[[0.0, 0.0, 0.0000], charge=-1, multiplicity=1) -basis_custom = get_basis_from_ccRepo(molecule, 'cc-pVTZ') +basis_custom = get_basis_from_BSE(molecule, 'cc-pVTZ') qc_input = QchemInput(molecule, jobtype='sp', diff --git a/examples/custom_basis.py b/examples/custom_basis.py index 2e38089..9bfcc3f 100644 --- a/examples/custom_basis.py +++ b/examples/custom_basis.py @@ -3,6 +3,7 @@ from pyqchem.structure import Structure from pyqchem.parsers.basic import basic_parser_qchem from pyqchem.basis import get_basis_from_ccRepo +from pyqchem.basis import get_basis_from_BSE # create molecule @@ -63,3 +64,19 @@ ) print('scf_energy (custom basis: cc-pVTZ): ', output['scf_energy']) + +basis_custom_repo = get_basis_from_BSE(molecule, '4zapa-nr', if_missing='6-31g') + +qc_input = QchemInput(molecule, + jobtype='sp', + exchange='hf', + basis=basis_custom_repo) + + +output = get_output_from_qchem(qc_input, + processors=4, + force_recalculation=False, + parser=basic_parser_qchem + ) + +print('scf_energy (custom basis: coemd): ', output['scf_energy']) diff --git a/pyqchem/__init__.py b/pyqchem/__init__.py index 42fd67b..f03d805 100644 --- a/pyqchem/__init__.py +++ b/pyqchem/__init__.py @@ -1,5 +1,5 @@ __author__ = 'Abel Carreras' -__version__ = '0.8.0' +__version__ = '0.8.1' from pyqchem.structure import Structure from pyqchem.qc_input import QchemInput diff --git a/pyqchem/basis.py b/pyqchem/basis.py index fd9d578..2ee6ef7 100644 --- a/pyqchem/basis.py +++ b/pyqchem/basis.py @@ -65,7 +65,11 @@ def get_basis_element_from_ccRepo(element, element_dict = {atom[1]: atom[2].lower() for atom in atom_data[1:]} - resp = req.get("http://www.grant-hill.group.shef.ac.uk/ccrepo/{}".format(element_dict[element])) + try: + resp = req.get("http://www.grant-hill.group.shef.ac.uk/ccrepo/{}".format(element_dict[element])) + except KeyError as e: + raise Exception('Atom label {} not recognized'.format(e.args[0])) + n_ini = resp.text.find('form-inline') # determinte the php web address for particular element @@ -93,6 +97,15 @@ def get_basis_element_from_ccRepo(element, def get_basis_from_ccRepo(structure, basis, full=False, if_missing=None): + """ + Get basis from ccRepo + + :param structure: Structure + :param basis: basis set label (string) + :param full: if False only list basis for unique atoms + :param if_missing: backup basis to use if basis is missing for a particular atom + :return: basis set dictionary + """ symbols = structure.get_symbols() if not full: @@ -110,9 +123,9 @@ def get_basis_from_ccRepo(structure, basis, full=False, if_missing=None): program='Gaussian', basis=if_missing) if len(basis_data) == 0: - raise Exception('Basis {}, {} not found for atom {}'.format(basis, if_missing, symbol)) + raise Exception('Basis {}, {} not found for atom {} in ccRepo'.format(basis, if_missing, symbol)) else: - raise Exception('Basis {} not found for atom {}'.format(basis, symbol)) + raise Exception('Basis {} not found for atom {} in ccRepo'.format(basis, symbol)) # print(description) # print(citation) @@ -125,6 +138,76 @@ def get_basis_from_ccRepo(structure, basis, full=False, if_missing=None): return basis_set +def get_basis_from_BSE(structure, basis, full=False, if_missing=None): + """ + get basis from Basis Set Exchange + + :param structure: + :param basis: basis set label (string) + :param full: if False only list basis for unique atoms + :return: basis set dictionary + """ + + base_url = "http://basissetexchange.org" + basis = basis.lower() + + headers = { + 'User-Agent': 'BSE Example Python Script', + 'From': 'bse@molssi.org' + } + + r = req.get(base_url + '/api/basis/{}/format/gaussian94'.format(basis), + headers=headers + ) + + description = '\n'.join([line for line in r.text.split('\n') if '!' in line]) + + r = req.get(base_url + '/api/references/{}/format/bib'.format(basis), + headers=headers + ) + + citation = r.text + + symbols = structure.get_symbols() + if not full: + symbols = np.unique(symbols) + + atoms = [] + for symbol in symbols: + + params = {'elements': [symbol]} + r = req.get(base_url + '/api/basis/{}/format/gaussian94'.format(basis), + params=params, + headers=headers + ) + # https://www.basissetexchange.org/basis/coemd-2/format/gaussian94/?version=0&elements=6,7,8 + if r.status_code != 200: + if if_missing is not None: + r = req.get(base_url + '/api/basis/{}/format/gaussian94'.format(if_missing), + params=params, + headers=headers + ) + if r.status_code != 200: + raise Exception('Basis {}, {} not found for atom {} in BSE'.format(basis, if_missing, symbol)) + else: + #raise RuntimeError("Could not obtain data from the BSE. Check the error information above") + raise Exception('Basis {} not found for atom {} in BSE'.format(basis, symbol)) + + basis_data = [] + for line in r.text.split('\n'): + if len(line) !=0 and line[0] not in ['!']: + basis_data.append(line.replace('D+', 'E+').replace('D-', 'E-')) + + + atoms.append(_txt_to_basis_dict(basis_data)) + + basis_set = {'name': basis, + 'primitive_type': 'gaussian', + 'atoms': atoms} + + return basis_set + + def basis_to_txt(basis): # write basis in qchem/gaussian format basis_txt = '' @@ -189,3 +272,12 @@ def trucate_basis(basis, shells=()): full=False) print(basis_to_txt(basis)) + + print('----------------------') + exit() + basis = get_basis_from_BSE(molecule, + basis='cc-pVTZ', + full=False) + + print(basis_to_txt(basis)) +