Secure your code as it's written. Use Snyk Code to scan source code in minutes - no build needed - and fix issues immediately.
def harvest_GRD(grd):
"""Parses the contents *grd* of the Cfour GRD file into the gradient
array and coordinate information. The coordinate info is converted
into a rather dinky Molecule (no charge, multiplicity, or fragment),
but this is these coordinates that govern the reading of molecule
orientation by Cfour. Return qcel.models.Molecule and gradient array.
"""
grd = grd.splitlines()
Nat = int(grd[0].split()[0])
molxyz = f"{Nat} bohr\n\n"
grad = []
for at in range(Nat):
mline = grd[at + 1].split()
el = "GH" if int(float(mline[0])) == 0 else qcel.periodictable.to_E(int(float(mline[0])))
molxyz += "%s %16s %16s %16s\n" % (el, mline[-3], mline[-2], mline[-1])
lline = grd[at + 1 + Nat].split()
grad.append([float(lline[-3]), float(lline[-2]), float(lline[-1])])
mol = Molecule(
validate=False,
**qcel.molparse.to_schema(
qcel.molparse.from_string(molxyz, dtype="xyz+", fix_com=True, fix_orientation=True)["qm"], dtype=2
),
)
return mol, grad
coords, number, _ = log.load_geometry()
except LogError:
logger.debug(f'Could not parse xyz from {path}')
# try parsing Gaussian standard orientation instead of the input orientation parsed by Arkane
lines = _get_lines_from_file(path)
xyz_str = ''
for i in range(len(lines)):
if 'Standard orientation:' in lines[i]:
xyz_str = ''
j = i
while not lines[j].split()[0].isdigit():
j += 1
while '-------------------' not in lines[j]:
splits = lines[j].split()
xyz_str += f'{qcel.periodictable.to_E(int(splits[1]))} {splits[3]} {splits[4]} {splits[5]}\n'
j += 1
break
if xyz_str:
return str_to_xyz(xyz_str)
return None
return xyz_from_data(coords=coords, numbers=number)
'Sc': 0.75, 'Ti': 0.86, 'V' : 0.79, 'Cr': 0.73, 'Mn': 0.67,
'Fe': 0.61, 'Co': 0.64, 'Ni': 0.55, 'Cu': 0.46, 'Zn': 0.60,
'Ga': 1.22, 'Ge': 1.22, 'As': 1.22, 'Se': 1.17, 'Br': 1.14, 'Kr': 1.03,
'I' : 1.33,
'X' : 0.00} # yapf: disable
#'RN': 2.40 / 1.5, # extrapolation
#'H': 1.06 / 1.5, # Bondi JPC 68 441 (1964)
#'SN': 2.16 / 1.5, # Bondi JPC 68 441 (1964)
#'SB': 2.12 / 1.5, # Bondi JPC 68 441 (1964)
#'TE': 2.08 / 1.5, # Bondi JPC 68 441 (1964)
#'XE': 2.05 / 1.5} # Bondi JPC 68 441 (1964)
nat = elem.shape[0]
try:
caps = [el.capitalize() for el in elem]
except AttributeError:
caps = [qcel.periodictable.to_E(z) for z in elem]
covrad = np.fromiter((covalent_radii_lookup[caps[at]] for at in range(nat)), dtype=np.float, count=nat)
return np.divide(covrad, qcel.constants.bohr2angstroms)
def sym2num(self, sym):
"""
Given a chemical symbol, returns the atomic number defined within the class
Parameters
-----------
sym : string
chemical symbol
Returns
--------
at_num : int
atomic number for symbol argument
"""
try:
atomic_number = qcel.periodictable.to_Z(sym)
return atomic_number
except:
raise KeyError('{} is not defined.'.format(sym))
def get_spinmult(self, symbol, charge=0):
anumber = qcelemental.periodictable.to_atomic_number(symbol)
nele = anumber - charge
spinmult = 0
#Neutral atoms
if charge == 0:
if anumber <= len(self.mult_dict):
spinmult = self.mult_dict[nele]
else:
raise ValueError("Determination of ion spinmult hasn't been\
implemented for elements with electrons occupying\
4d orbitals")
else:
if anumber <= 20: # H to Ca
spinmult = self.mult_dict[nele]
elif anumber <= 30: #Sc - Zn
if charge == 1:#lose one 4s electron
# Assume that electrons are always removed from 4s
def offer_mass_number(z, a):
"""Given a mass number and element, what can be suggested and asserted about A, mass?"""
a = int(a)
a_eliso = qcel.periodictable.to_E(z) + str(a)
a_mass = qcel.periodictable.to_mass(a_eliso)
A_exact.append(a)
A_range.append(lambda x, a=a: x == a)
text.append("""For A, inputs Z: {}, A: {} require A == {}.""".format(z, a, a))
m_exact.append(a_mass)
m_range.append(lambda x, a_mass=a_mass: abs(x - a_mass) < mtol)
text.append("""For mass, inputs Z: {}, A: {} require abs(mass - {}) < {}""".format(z, a, a_mass, mtol))
molrec['fix_symmetry'] = 'c1'
elif self.symmetry_from_input():
molrec['fix_symmetry'] = self.symmetry_from_input()
# if self.has_zmatrix:
# moldict['zmat'] = self.zmat
# TODO zmat, geometry_variables
nat = self.natom()
geom = np.array(self.geometry()) # [a0]
if molrec['units'] == 'Angstrom':
geom *= qcel.constants.bohr2angstroms #self.input_units_to_au()
molrec['geom'] = geom.reshape((-1))
molrec['elea'] = np.array([self.mass_number(at) for at in range(nat)])
molrec['elez'] = np.array([qcel.periodictable.to_Z(self.symbol(at)) for at in range(nat)])
molrec['elem'] = np.array([self.symbol(at).capitalize() for at in range(nat)])
molrec['mass'] = np.array([self.mass(at) for at in range(nat)])
molrec['real'] = np.array([bool(self.Z(at)) for at in range(nat)])
molrec['elbl'] = np.array([self.label(at)[len(self.symbol(at)):].lower() for at in range(nat)])
fragments = [x[:] for x in self.get_fragments()]
fragment_charges = [float(f) for f in self.get_fragment_charges()]
fragment_multiplicities = [m for m in self.get_fragment_multiplicities()]
# do trimming not performed in Molecule class b/c fragment_* member data never directly exposed
for ifr, fr in reversed(list(enumerate(self.get_fragment_types()))):
if fr == 'Ghost':
fragment_charges[ifr] = 0.
fragment_multiplicities[ifr] = 1
elif fr == 'Absent':
del fragment_charges[ifr]
def mass(self, atom):
"""Returns mass of atom (0-indexed)
>>> print(H2OH2O.mass(4))
1.00782503207
"""
if self.atoms[atom].mass() != 0.0:
return self.atoms[atom].mass()
if math.fabs(self.atoms[atom].Z() - int(self.atoms[atom].Z())) > 0.0:
print("""WARNING: Obtaining masses from atom with fractional charge...may be incorrect!!!\n""")
# TODO outfile
return qcel.periodictable.to_mass(int(self.atoms[atom].Z()))
def send_elements(self):
""" Send the atomic number of each nucleus through MDI
Returns
-------
elements : :obj:`list` of :obj:`int`
Element of each atom
"""
natom = len(self.molecule.geometry)
elements = [qcel.periodictable.to_atomic_number(self.molecule.symbols[iatom]) for iatom in range(natom)]
MDI_Send(elements, natom, MDI_INT, self.comm)
return elements