Skip to content

Commit

Permalink
Merge remote-tracking branch 'origin/master'
Browse files Browse the repository at this point in the history
  • Loading branch information
@abelcarreras
abelcarreras committed Mar 10, 2021
2 parents 18ee409 + 36bc210 commit 7a884cc
Show file tree
Hide file tree
Showing 5 changed files with 142 additions and 13 deletions.
14 changes: 9 additions & 5 deletions cosymlib/__init__.py
View file
Original file line number Diff line number Diff line change
@@ -1,4 +1,4 @@
__version__ = '0.9.4'
__version__ = '0.9.5'

from cosymlib.molecule import Molecule
from cosymlib.molecule.geometry import Geometry
Expand Down Expand Up @@ -509,7 +509,8 @@ def print_wnfsym_sym_ovelap(self, group, axis=None, axis2=None, center=None, out

output.write(txt)

def print_wnfsym_irreducible_repr(self, group, axis=None, axis2=None, center=None, output=sys.stdout):
def print_wnfsym_irreducible_repr(self, group, axis=None, axis2=None, center=None, show_axes=True,
output=sys.stdout):

txt = ''
for molecule in self._molecules:
Expand All @@ -531,11 +532,13 @@ def print_wnfsym_irreducible_repr(self, group, axis=None, axis2=None, center=Non
txt += 'WFN ' + ' '.join(['{:7.3f}'.format(s) for s in data_wf['total']])
txt += '\n'

txt += _get_axis_info(molecule, group, axis, axis2, center)
if show_axes:
txt += _get_axis_info(molecule, group, axis, axis2, center)

output.write(txt)

def print_wnfsym_mo_irreducible_repr(self, group, axis=None, axis2=None, center=None, output=sys.stdout):
def print_wnfsym_mo_irreducible_repr(self, group, axis=None, axis2=None, center=None, show_axes=True,
output=sys.stdout):

txt = ''
for molecule in self._molecules:
Expand Down Expand Up @@ -581,7 +584,8 @@ def print_wnfsym_mo_irreducible_repr(self, group, axis=None, axis2=None, center=
if non_one:
warn('The sum of some molecular orbital irreducible represaentations are not equal one (*)')

txt += _get_axis_info(molecule, group, axis, axis2, center)
if show_axes:
txt += _get_axis_info(molecule, group, axis, axis2, center)

output.write(txt)

Expand Down
133 changes: 131 additions & 2 deletions cosymlib/file_io/__init__.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -2,7 +2,7 @@
from cosymlib.molecule.geometry import Geometry
from cosymlib.molecule.electronic_structure import ElectronicStructure
from cosymlib.file_io.tools import extract_geometries
from cosymlib.tools import atomic_number_to_element
from cosymlib.tools import atomic_number_to_element, element_to_atomic_number
from cosymlib.file_io import errors

import numpy as np
Expand Down Expand Up @@ -72,7 +72,7 @@ def get_geometry_from_file_xyz(file_name, read_multiple=False):
float(line.split()[1])
input_molecule[0].append(line.split()[0])
input_molecule[1].append(line.split()[1:])
except IndexError:
except (IndexError, ValueError):
if input_molecule[0]:
if len(input_molecule[0]) != n_atoms:
warnings.warn('Number of atoms in line {} and number '
Expand Down Expand Up @@ -574,3 +574,132 @@ def reformat_input(array):
else:
flat_list.append(item)
return flat_list


def get_array_txt(label, type, array, row_size=5):

formats = {'R': '15.8e',
'I': '11'}

n_elements = len(array)
rows = int(np.ceil(n_elements/row_size))

txt_fchk = '{:40} {} N= {:5}\n'.format(label, type, n_elements)

for i in range(rows):
if (i+1)*row_size > n_elements:
txt_fchk += (' {:{fmt}}'* (n_elements - i*row_size) + '\n').format(*array[i * row_size:n_elements],
fmt=formats[type])
else:
txt_fchk += (' {:{fmt}}'* row_size + '\n').format(*array[i * row_size: (i+1)*row_size],
fmt=formats[type])

return txt_fchk


def build_fchk(molecule):


structure = np.array(molecule.geometry.get_positions())
basis = molecule.electronic_structure.basis
alpha_mo_coeff = molecule.electronic_structure.coefficients_a
alpha_mo_energies = molecule.electronic_structure.alpha_energies
if molecule.electronic_structure.coefficients_b != molecule.electronic_structure.coefficients_a:
beta_mo_coeff = molecule.electronic_structure.coefficients_b
beta_mo_energies = molecule.electronic_structure.beta_energies
beta_mo_coeff = np.array(beta_mo_coeff).flatten().tolist()
else:
beta_mo_coeff = None
beta_mo_energies = None

number_of_basis_functions = len(alpha_mo_coeff)
alpha_electrons = molecule.electronic_structure.alpha_electrons
beta_electrons = molecule.electronic_structure.beta_electrons
number_of_electrons = alpha_electrons + beta_electrons

alpha_mo_coeff = np.array(alpha_mo_coeff).flatten().tolist()

atomic_numbers = [element_to_atomic_number(symbol) for symbol in molecule.geometry.get_symbols()]
total_energy = np.dot(alpha_mo_energies, alpha_electrons) + np.dot(beta_mo_energies, beta_electrons)

shell_type_list = {'s': {'type': 0, 'angular_momentum': 0},
'p': {'type': 1, 'angular_momentum': 1},
'd': {'type': 2, 'angular_momentum': 2},
'f': {'type': 3, 'angular_momentum': 3},
'sp': {'type': -1, 'angular_momentum': 1}, # hybrid
'd_': {'type': -2, 'angular_momentum': 2}, # pure
'f_': {'type': -3, 'angular_momentum': 3}} # pure

shell_type = []
p_exponents = []
c_coefficients = []
p_c_coefficients = []
n_primitives = []
atom_map = []

largest_degree_of_contraction = 0
highest_angular_momentum = 0
number_of_contracted_shells = 0

for i, atoms in enumerate(basis['atoms']):
for shell in atoms['shells']:
number_of_contracted_shells += 1
st = shell['shell_type']
shell_type.append(shell_type_list[st]['type'])
n_primitives.append(len(shell['p_exponents']))
atom_map.append(i+1)
if highest_angular_momentum < shell_type_list[st]['angular_momentum']:
highest_angular_momentum = shell_type_list[st]['angular_momentum']

if len(shell['con_coefficients']) > largest_degree_of_contraction:
largest_degree_of_contraction = len(shell['con_coefficients'])

for p in shell['p_exponents']:
p_exponents.append(p)
for c in shell['con_coefficients']:
c_coefficients.append(c)
for pc in shell['p_con_coefficients']:
p_c_coefficients.append(pc)

angstrom_to_bohr = 1/0.529177249
coordinates_list = angstrom_to_bohr*structure.flatten()

txt_fchk = '{}\n'.format('Extended Huckel calculation over filename')
txt_fchk += 'SP RHuckel {}\n'.format('STO-3G')
txt_fchk += 'Number of atoms I {}\n'.format(len(structure))
txt_fchk += 'Charge I {}\n'.format(molecule.get_charge())
txt_fchk += 'Multiplicity I {}\n'.format(molecule.electronic_structure.
multiplicity)
txt_fchk += 'Number of electrons I {}\n'.format(number_of_electrons)
txt_fchk += 'Number of alpha electrons I {}\n'.format(alpha_electrons)
txt_fchk += 'Number of beta electrons I {}\n'.format(beta_electrons)
txt_fchk += 'Number of basis functions I {}\n'.format(number_of_basis_functions)

txt_fchk += get_array_txt('Atomic numbers', 'I', atomic_numbers, row_size=6)
txt_fchk += get_array_txt('Nuclear charges', 'R', atomic_numbers)
txt_fchk += get_array_txt('Current cartesian coordinates', 'R', coordinates_list)

txt_fchk += 'Number of contracted shells I {}\n'.format(number_of_contracted_shells)
txt_fchk += 'Number of primitive shells I {}\n'.format(np.sum(n_primitives))
txt_fchk += 'Highest angular momentum I {}\n'.format(highest_angular_momentum)
txt_fchk += 'Largest degree of contraction I {}\n'.format(largest_degree_of_contraction)

txt_fchk += get_array_txt('Shell types', 'I', shell_type, row_size=6)
txt_fchk += get_array_txt('Number of primitives per shell', 'I', n_primitives, row_size=6)
txt_fchk += get_array_txt('Shell to atom map', 'I', atom_map, row_size=6)
txt_fchk += get_array_txt('Primitive exponents', 'R', p_exponents)
txt_fchk += get_array_txt('Contraction coefficients', 'R', c_coefficients)
txt_fchk += get_array_txt('P(S=P) Contraction coefficients', 'R', p_c_coefficients)
# txt_fchk += get_array_txt('Coordinates of each shell', 'R', coor_shell) #
# txt_fchk += get_array_txt('Overlap Matrix', 'R', overlap)
#txt_fchk += get_array_txt('Core Hamiltonian Matrix', 'R', core_hamiltonian)
txt_fchk += 'Total Energy R {}\n'.format(total_energy)
txt_fchk += get_array_txt('Alpha Orbital Energies', 'R', alpha_mo_energies)
if beta_mo_energies is not None:
txt_fchk += get_array_txt('Beta Orbital Energies', 'R', beta_mo_energies)
# txt_fchk += get_array_txt('Total SCF Density', 'R', scf_density)
txt_fchk += get_array_txt('Alpha MO coefficients', 'R', alpha_mo_coeff)
if beta_mo_coeff is not None:
txt_fchk += get_array_txt('Beta MO coefficients', 'R', beta_mo_coeff)

return txt_fchk
6 changes: 0 additions & 6 deletions cosymlib/molecule/electronic_structure/__init__.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -26,14 +26,12 @@ class ElectronicStructure:
def __init__(self,
basis,
orbital_coefficients,
# charge=None,
multiplicity=None,
alpha_energies=None,
beta_energies=None,
alpha_occupancy=None,
beta_occupancy=None):

# self._charge = charge
self._multiplicity = multiplicity
self._basis = basis
self._Ca = orbital_coefficients[0]
Expand Down Expand Up @@ -88,10 +86,6 @@ def _occupancy_consistency(self):
for i in range(len(self._Cb) - len(self._beta_occupancy)):
self._beta_occupancy.append(0)

# @property
# def charge(self):
# return self._charge

@property
def multiplicity(self):
return self._multiplicity
Expand Down
1 change: 1 addition & 0 deletions scripts/cosym
View file
Original file line number Diff line number Diff line change
Expand Up @@ -265,6 +265,7 @@ if args.qsym_wfn:
axis=args.axis,
axis2=args.axis2,
center=args.center,
show_axes=False,
output=common_output)

symobj.print_wnfsym_irreducible_repr(args.qsym_wfn,
Expand Down
1 change: 1 addition & 0 deletions scripts/esym
View file
Original file line number Diff line number Diff line change
Expand Up @@ -147,6 +147,7 @@ if args.measure:
axis=args.axis,
axis2=args.axis2,
center=args.center,
show_axes=False,
output=common_output)

structure_set.print_wnfsym_irreducible_repr(args.measure,
Expand Down

0 comments on commit 7a884cc

Please sign in to comment.