Secure your code as it's written. Use Snyk Code to scan source code in minutes - no build needed - and fix issues immediately.
-0.0153816024 -0.0172298679 0.0000000000
0.1484742948 0.3028992171 0.0000000000
0.0000000000 0.0000000000 0.0073781116
0.0000000000 0.0000000000 -0.0450044998
0.0000000000 0.0000000000 0.0146744312
0.0000000000 0.0000000000 -0.0099324434
0.0000000000 0.0000000000 0.0041387542
0.0000000000 0.0000000000 0.0287456462
"""
ref_eth_vibonly = {
# freq from cfour after proj
'omega': qcel.Datum('', '', np.asarray([885.4386, 1073.4306, 1080.4104, 1135.9390, 1328.7166, 1467.8428, 1565.5328, 1831.6445, 3287.3173, 3312.1889, 3371.7974, 3399.4232])),
'gamma': qcel.Datum('', '', ['B2u', 'B1u', 'B2g', 'Au', 'B1g', 'Ag', 'B3u', 'Ag', 'B3u', 'Ag', 'B1g', 'B2u'], numeric=False),
'IR_intensity': qcel.Datum('', '', np.array([0.2122, 94.4517, 0.0000, 0.0000, 0.0000, 0.0000, 10.3819, 0.0000, 19.2423, 0.0000, 0.0000, 28.0922])), # [km/mol]
# from vibsuite after cfour
'mu': qcel.Datum('', '', np.asarray([1.0423, 1.1607, 1.5200, 1.0078, 1.5276, 1.2029, 1.1122, 3.2209, 1.0475, 1.0799, 1.1150, 1.1182])),
'k': qcel.Datum('', '', np.asarray([0.4814, 0.7880, 1.0454, 0.7662, 1.5890, 1.5270, 1.6060, 6.3666, 6.6695, 6.9800, 7.4688, 7.6133])),
'Qtp0': qcel.Datum('', '', np.asarray([0.3687, 0.3349, 0.3338, 0.3256, 0.3010, 0.2864, 0.2773, 0.2564, 0.1914, 0.1907, 0.1890, 0.1882])),
'DQ0': qcel.Datum('', '', np.asarray([0.2607, 0.2368, 0.2360, 0.2302, 0.2128, 0.2025, 0.1961, 0.1813, 0.1353, 0.1348, 0.1336, 0.1331])),
# q from vibsuite after cfour
'q': qcel.Datum('', '', np.asarray(
[[ 0.00000000E+00, -1.34306554E-01, 0.00000000E+00, 0.00000000E+00, -1.34306554E-01, 0.00000000E+00, -4.32765855E-01, 2.31720930E-01, 0.00000000E+00, 4.32765855E-01, 2.31720930E-01, 0.00000000E+00, 4.32765855E-01, 2.31720930E-01, 0.00000000E+00, -4.32765855E-01, 2.31720930E-01, 0.00000000E+00],
[ 0.00000000E+00, 0.00000000E+00, -2.68155030E-01, 0.00000000E+00, 0.00000000E+00, -2.68155030E-01, 0.00000000E+00, 0.00000000E+00, 4.62651532E-01, 0.00000000E+00, 0.00000000E+00, 4.62651532E-01, 0.00000000E+00, 0.00000000E+00, 4.62651532E-01, 0.00000000E+00, 0.00000000E+00, 4.62651532E-01],
[ 0.00000000E+00, 0.00000000E+00, -4.28876240E-01, 0.00000000E+00, 0.00000000E+00, 4.28876240E-01, 0.00000000E+00, 0.00000000E+00, 3.97533125E-01, 0.00000000E+00, 0.00000000E+00, 3.97533125E-01, 0.00000000E+00, 0.00000000E+00, -3.97533125E-01, 0.00000000E+00, 0.00000000E+00, -3.97533125E-01],
[ 0.00000000E+00, 0.00000000E+00, 0.00000000E+00, 0.00000000E+00, 0.00000000E+00, 0.00000000E+00, 0.00000000E+00, 0.00000000E+00, -5.00000000E-01, 0.00000000E+00, 0.00000000E+00, 5.00000000E-01, 0.00000000E+00, 0.00000000E+00, 5.00000000E-01, 0.00000000E+00, 0.00000000E+00, -5.00000000E-01],
[ 0.00000000E+00, 4.30954858E-01, 0.00000000E+00, 0.00000000E+00, -4.30954858E-01, 0.00000000E+00, 3.79464410E-01, -1.14654771E-01, 0.00000000E+00, -3.79464410E-01, -1.14654771E-01, 0.00000000E+00, 3.79464410E-01, 1.14654771E-01, 0.00000000E+00, -3.79464410E-01, 1.14654771E-01, 0.00000000E+00],
[ -2.97546697E-01, 0.00000000E+00, 0.00000000E+00, 2.97546697E-01, 0.00000000E+00, 0.00000000E+00, -4.15835441E-01, 1.81145984E-01, 0.00000000E+00, -4.15835441E-01, -1.81145984E-01, 0.00000000E+00, 4.15835441E-01, 1.81145984E-01, 0.00000000E+00, 4.15835441E-01, -1.81145984E-01, 0.00000000E+00],
[ 2.26323959E-01, 0.00000000E+00, 0.00000000E+00, 2.26323959E-01, 0.00000000E+00, 0.00000000E+00, -3.90479813E-01, 2.68168321E-01, 0.00000000E+00, -3.90479813E-01, -2.68168321E-01, 0.00000000E+00, -3.90479813E-01, -2.68168321E-01, 0.00000000E+00, -3.90479813E-01, 2.68168321E-01, 0.00000000E+00],
'E_elec' : qcel.Datum('', '', 0.00000000),
'E_trans' : qcel.Datum('', '', 0.00141628),
'E_rot' : qcel.Datum('', '', 0.00141628),
'E_vib' : qcel.Datum('', '', 0.04801283),
'E_corr' : qcel.Datum('', '', 0.05084539),
'E_tot' : qcel.Datum('', '', -39.92603216),
'H_elec' : qcel.Datum('', '', 0.00000000),
'H_trans' : qcel.Datum('', '', 0.00236046),
'H_rot' : qcel.Datum('', '', 0.00141628),
'H_vib' : qcel.Datum('', '', 0.04801283),
'H_corr' : qcel.Datum('', '', 0.05178958),
'H_tot' : qcel.Datum('', '', -39.92508798),
'G_elec' : qcel.Datum('', '', -0.00000000),
'G_trans' : qcel.Datum('', '', -0.01391819),
'G_rot' : qcel.Datum('', '', -0.00337341),
'G_vib' : qcel.Datum('', '', 0.04799343),
'G_corr' : qcel.Datum('', '', 0.03070183),
'G_tot' : qcel.Datum('', '', -39.94617572),
} # yapf: disable
assert qcdb.compare_vibinfos(ch4_hf_321g_thermoinfo, therminfo, 4, 'asdf', forgive=['omega', 'IR_intensity'])
try:
fullgrad[ireal, :] = rg
except NameError as err:
raise Dftd3Error('Unsuccessful gradient collection.') from err
qcvkey = jobrec['options']['fctldash'].upper()
# OLD WAY
calcinfo = []
if dftd3rec['dashlevel'] == 'atmgr':
calcinfo.append(qcel.Datum('DISPERSION CORRECTION ENERGY', 'Eh', atm))
calcinfo.append(qcel.Datum('3-BODY DISPERSION CORRECTION ENERGY', 'Eh', atm))
calcinfo.append(qcel.Datum('AXILROD-TELLER-MUTO 3-BODY DISPERSION CORRECTION ENERGY', 'Eh', atm))
if jobrec['driver'] == 'gradient':
calcinfo.append(qcel.Datum('DISPERSION CORRECTION GRADIENT', 'Eh/a0', fullgrad))
calcinfo.append(qcel.Datum('3-BODY DISPERSION CORRECTION GRADIENT', 'Eh/a0', fullgrad))
calcinfo.append(qcel.Datum('AXILROD-TELLER-MUTO 3-BODY DISPERSION CORRECTION GRADIENT', 'Eh/a0', fullgrad))
else:
calcinfo.append(qcel.Datum('DISPERSION CORRECTION ENERGY', 'Eh', ene))
calcinfo.append(qcel.Datum('2-BODY DISPERSION CORRECTION ENERGY', 'Eh', ene))
if qcvkey:
calcinfo.append(qcel.Datum(f'{qcvkey} DISPERSION CORRECTION ENERGY', 'Eh', ene))
if jobrec['driver'] == 'gradient':
calcinfo.append(qcel.Datum('DISPERSION CORRECTION GRADIENT', 'Eh/a0', fullgrad))
calcinfo.append(qcel.Datum('2-BODY DISPERSION CORRECTION GRADIENT', 'Eh/a0', fullgrad))
if qcvkey:
calcinfo.append(qcel.Datum(f'{qcvkey} DISPERSION CORRECTION GRADIENT', 'Eh/a0', fullgrad))
calcinfo = {info.label: info for info in calcinfo}
qcvkey = input_model.extras["info"]["fctldash"].upper()
calcinfo = []
calcinfo.append(qcel.Datum("CURRENT ENERGY", "Eh", ene))
calcinfo.append(qcel.Datum("DISPERSION CORRECTION ENERGY", "Eh", ene))
calcinfo.append(qcel.Datum("2-BODY DISPERSION CORRECTION ENERGY", "Eh", ene))
if qcvkey:
calcinfo.append(qcel.Datum(f"{qcvkey} DISPERSION CORRECTION ENERGY", "Eh", ene))
if input_model.driver == "gradient":
calcinfo.append(qcel.Datum("CURRENT GRADIENT", "Eh/a0", fullgrad))
calcinfo.append(qcel.Datum("DISPERSION CORRECTION GRADIENT", "Eh/a0", fullgrad))
calcinfo.append(qcel.Datum("2-BODY DISPERSION CORRECTION GRADIENT", "Eh/a0", fullgrad))
if qcvkey:
calcinfo.append(qcel.Datum(f"{qcvkey} DISPERSION CORRECTION GRADIENT", "Eh/a0", fullgrad))
# LOGtext += qcel.datum.print_variables({info.label: info for info in calcinfo})
calcinfo = {info.label: info.data for info in calcinfo}
# calcinfo = qcel.util.unnp(calcinfo, flat=True)
# got to even out who needs plump/flat/Decimal/float/ndarray/list
# Decimal --> str preserves precision
calcinfo = {
k.upper(): str(v) if isinstance(v, Decimal) else v for k, v in qcel.util.unnp(calcinfo, flat=True).items()
}
# jobrec['properties'] = {"return_energy": ene}
# jobrec["molecule"]["real"] = list(jobrec["molecule"]["real"])
retres = calcinfo[f"CURRENT {input_model.driver.upper()}"]
if isinstance(retres, Decimal):
try:
fullgrad[ireal, :] = realgrad
except NameError as exc:
raise UnknownError("Unsuccessful gradient collection.") from exc
qcvkey = input_model.extras["info"]["fctldash"].upper()
calcinfo = []
calcinfo.append(qcel.Datum("CURRENT ENERGY", "Eh", ene))
calcinfo.append(qcel.Datum("DISPERSION CORRECTION ENERGY", "Eh", ene))
calcinfo.append(qcel.Datum("2-BODY DISPERSION CORRECTION ENERGY", "Eh", ene))
if qcvkey:
calcinfo.append(qcel.Datum(f"{qcvkey} DISPERSION CORRECTION ENERGY", "Eh", ene))
if input_model.driver == "gradient":
calcinfo.append(qcel.Datum("CURRENT GRADIENT", "Eh/a0", fullgrad))
calcinfo.append(qcel.Datum("DISPERSION CORRECTION GRADIENT", "Eh/a0", fullgrad))
calcinfo.append(qcel.Datum("2-BODY DISPERSION CORRECTION GRADIENT", "Eh/a0", fullgrad))
if qcvkey:
calcinfo.append(qcel.Datum(f"{qcvkey} DISPERSION CORRECTION GRADIENT", "Eh/a0", fullgrad))
# LOGtext += qcel.datum.print_variables({info.label: info for info in calcinfo})
calcinfo = {info.label: info.data for info in calcinfo}
# calcinfo = qcel.util.unnp(calcinfo, flat=True)
# got to even out who needs plump/flat/Decimal/float/ndarray/list
# Decimal --> str preserves precision
calcinfo = {
k.upper(): str(v) if isinstance(v, Decimal) else v for k, v in qcel.util.unnp(calcinfo, flat=True).items()
}
# jobrec['properties'] = {"return_energy": ene}
calcinfo.append(qcel.Datum('AXILROD-TELLER-MUTO 3-BODY DISPERSION CORRECTION ENERGY', 'Eh', atm))
if jobrec['driver'] == 'gradient':
calcinfo.append(qcel.Datum('DISPERSION CORRECTION GRADIENT', 'Eh/a0', fullgrad))
calcinfo.append(qcel.Datum('3-BODY DISPERSION CORRECTION GRADIENT', 'Eh/a0', fullgrad))
calcinfo.append(qcel.Datum('AXILROD-TELLER-MUTO 3-BODY DISPERSION CORRECTION GRADIENT', 'Eh/a0', fullgrad))
else:
calcinfo.append(qcel.Datum('DISPERSION CORRECTION ENERGY', 'Eh', ene))
calcinfo.append(qcel.Datum('2-BODY DISPERSION CORRECTION ENERGY', 'Eh', ene))
if qcvkey:
calcinfo.append(qcel.Datum(f'{qcvkey} DISPERSION CORRECTION ENERGY', 'Eh', ene))
if jobrec['driver'] == 'gradient':
calcinfo.append(qcel.Datum('DISPERSION CORRECTION GRADIENT', 'Eh/a0', fullgrad))
calcinfo.append(qcel.Datum('2-BODY DISPERSION CORRECTION GRADIENT', 'Eh/a0', fullgrad))
if qcvkey:
calcinfo.append(qcel.Datum(f'{qcvkey} DISPERSION CORRECTION GRADIENT', 'Eh/a0', fullgrad))
calcinfo = {info.label: info for info in calcinfo}
text += qcel.datum.print_variables(calcinfo)
# NEW WAY
#module_vars = {}
#module_vars['DISPERSION CORRECTION ENERGY'] = ene
#module_vars['{} DISPERSION CORRECTION ENERGY'.format(qcvkey)] = ene
#if jobrec['driver'] == 'gradient':
# module_vars['DISPERSION CORRECTION GRADIENT'] = fullgrad
# module_vars['{} DISPERSION CORRECTION GRADIENT'.format(qcvkey)] = fullgrad
#
#module_vars = PreservingDict(module_vars)
#qcvars.build_out(module_vars)
mwhess_proj = np.dot(P.T, mwhess).dot(P)
text.append(mat_symm_info(mwhess_proj, lbl='projected mass-weighted Hessian') + f' ({nrt})')
#print('projhess = ', np.array_repr(mwhess_proj))
force_constant_au, qL = np.linalg.eigh(mwhess_proj)
# expected order for vibrations is steepest downhill to steepest uphill
idx = np.argsort(force_constant_au)
force_constant_au = force_constant_au[idx]
qL = qL[:, idx]
qL = _phase_cols_to_max_element(qL)
vibinfo['q'] = Datum('normal mode', 'a0 u^1/2', qL, comment='normalized mass-weighted')
# frequency, LAB II.17
frequency_cm_1 = np.lib.scimath.sqrt(force_constant_au) * uconv_cm_1
vibinfo['omega'] = Datum('frequency', 'cm^-1', frequency_cm_1)
# degeneracies
ufreq, uinv, ucts = np.unique(np.around(frequency_cm_1, 2), return_inverse=True, return_counts=True)
vibinfo['degeneracy'] = Datum('degeneracy', '', ucts[uinv])
# look among the symmetry subspaces h for one to which the normco
# of vib does *not* add an extra dof to the vector space
active = []
irrep_classification = []
for idx, vib in enumerate(frequency_cm_1):
if vec_in_space(qL[:, idx], TRspace, 1.0e-4):
active.append('TR')
irrep_classification.append(None)
else:
realgradabc = np.zeros((1, 3))
if input_model.driver == "gradient":
ireal = np.argwhere(real).reshape((-1))
fullgrad = np.zeros((full_nat, 3))
rg = realgradabc if (input_model.extras["info"]["dashlevel"] == "atmgr") else realgrad
try:
fullgrad[ireal, :] = rg
except NameError as exc:
raise UnknownError("Unsuccessful gradient collection.") from exc
qcvkey = input_model.extras["info"]["fctldash"].upper()
calcinfo = []
if input_model.extras["info"]["dashlevel"] == "atmgr":
calcinfo.append(qcel.Datum("CURRENT ENERGY", "Eh", atm))
calcinfo.append(qcel.Datum("DISPERSION CORRECTION ENERGY", "Eh", atm))
calcinfo.append(qcel.Datum("3-BODY DISPERSION CORRECTION ENERGY", "Eh", atm))
calcinfo.append(qcel.Datum("AXILROD-TELLER-MUTO 3-BODY DISPERSION CORRECTION ENERGY", "Eh", atm))
if input_model.driver == "gradient":
calcinfo.append(qcel.Datum("CURRENT GRADIENT", "Eh/a0", fullgrad))
calcinfo.append(qcel.Datum("DISPERSION CORRECTION GRADIENT", "Eh/a0", fullgrad))
calcinfo.append(qcel.Datum("3-BODY DISPERSION CORRECTION GRADIENT", "Eh/a0", fullgrad))
calcinfo.append(
qcel.Datum("AXILROD-TELLER-MUTO 3-BODY DISPERSION CORRECTION GRADIENT", "Eh/a0", fullgrad)
)
else:
calcinfo.append(qcel.Datum("CURRENT ENERGY", "Eh", ene))
calcinfo.append(qcel.Datum("DISPERSION CORRECTION ENERGY", "Eh", ene))
calcinfo.append(qcel.Datum("2-BODY DISPERSION CORRECTION ENERGY", "Eh", ene))
qcvkey = input_model.extras["info"]["fctldash"].upper()
calcinfo = []
if input_model.extras["info"]["dashlevel"] == "atmgr":
calcinfo.append(qcel.Datum("CURRENT ENERGY", "Eh", atm))
calcinfo.append(qcel.Datum("DISPERSION CORRECTION ENERGY", "Eh", atm))
calcinfo.append(qcel.Datum("3-BODY DISPERSION CORRECTION ENERGY", "Eh", atm))
calcinfo.append(qcel.Datum("AXILROD-TELLER-MUTO 3-BODY DISPERSION CORRECTION ENERGY", "Eh", atm))
if input_model.driver == "gradient":
calcinfo.append(qcel.Datum("CURRENT GRADIENT", "Eh/a0", fullgrad))
calcinfo.append(qcel.Datum("DISPERSION CORRECTION GRADIENT", "Eh/a0", fullgrad))
calcinfo.append(qcel.Datum("3-BODY DISPERSION CORRECTION GRADIENT", "Eh/a0", fullgrad))
calcinfo.append(
qcel.Datum("AXILROD-TELLER-MUTO 3-BODY DISPERSION CORRECTION GRADIENT", "Eh/a0", fullgrad)
)
else:
calcinfo.append(qcel.Datum("CURRENT ENERGY", "Eh", ene))
calcinfo.append(qcel.Datum("DISPERSION CORRECTION ENERGY", "Eh", ene))
calcinfo.append(qcel.Datum("2-BODY DISPERSION CORRECTION ENERGY", "Eh", ene))
if qcvkey:
calcinfo.append(qcel.Datum(f"{qcvkey} DISPERSION CORRECTION ENERGY", "Eh", ene))
if input_model.driver == "gradient":
calcinfo.append(qcel.Datum("CURRENT GRADIENT", "Eh/a0", fullgrad))
calcinfo.append(qcel.Datum("DISPERSION CORRECTION GRADIENT", "Eh/a0", fullgrad))
calcinfo.append(qcel.Datum("2-BODY DISPERSION CORRECTION GRADIENT", "Eh/a0", fullgrad))
if qcvkey:
calcinfo.append(qcel.Datum(f"{qcvkey} DISPERSION CORRECTION GRADIENT", "Eh/a0", fullgrad))
srealgrad = outfiles["mp2d_gradient"]
realgrad = np.fromstring(srealgrad, count=3 * real_nat, sep=" ").reshape((-1, 3))
if input_model.driver == "gradient":
ireal = np.argwhere(real).reshape((-1))
fullgrad = np.zeros((full_nat, 3))
try:
fullgrad[ireal, :] = realgrad
except NameError as exc:
raise UnknownError("Unsuccessful gradient collection.") from exc
qcvkey = input_model.extras["info"]["fctldash"].upper()
calcinfo = []
calcinfo.append(qcel.Datum("CURRENT ENERGY", "Eh", ene))
calcinfo.append(qcel.Datum("DISPERSION CORRECTION ENERGY", "Eh", ene))
calcinfo.append(qcel.Datum("2-BODY DISPERSION CORRECTION ENERGY", "Eh", ene))
if qcvkey:
calcinfo.append(qcel.Datum(f"{qcvkey} DISPERSION CORRECTION ENERGY", "Eh", ene))
if input_model.driver == "gradient":
calcinfo.append(qcel.Datum("CURRENT GRADIENT", "Eh/a0", fullgrad))
calcinfo.append(qcel.Datum("DISPERSION CORRECTION GRADIENT", "Eh/a0", fullgrad))
calcinfo.append(qcel.Datum("2-BODY DISPERSION CORRECTION GRADIENT", "Eh/a0", fullgrad))
if qcvkey:
calcinfo.append(qcel.Datum(f"{qcvkey} DISPERSION CORRECTION GRADIENT", "Eh/a0", fullgrad))
# LOGtext += qcel.datum.print_variables({info.label: info for info in calcinfo})
calcinfo = {info.label: info.data for info in calcinfo}
# calcinfo = qcel.util.unnp(calcinfo, flat=True)
# got to even out who needs plump/flat/Decimal/float/ndarray/list