Skip to content

Commit

Permalink
Added function for getting custom basis from BSE
Browse files Browse the repository at this point in the history
  • Loading branch information
abelcarreras committed Aug 13, 2021
1 parent 46b3346 commit 37c9df3
Show file tree
Hide file tree
Showing 4 changed files with 115 additions and 6 deletions.
4 changes: 2 additions & 2 deletions README.md
Original file line number Diff line number Diff line change
Expand Up @@ -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],
Expand All @@ -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',
Expand Down
17 changes: 17 additions & 0 deletions examples/custom_basis.py
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down Expand Up @@ -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'])
2 changes: 1 addition & 1 deletion pyqchem/__init__.py
Original file line number Diff line number Diff line change
@@ -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
Expand Down
98 changes: 95 additions & 3 deletions pyqchem/basis.py
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down Expand Up @@ -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:
Expand All @@ -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)
Expand All @@ -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 = ''
Expand Down Expand Up @@ -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))

0 comments on commit 37c9df3

Please sign in to comment.