Source code for diatomic

#
#@BEGIN LICENSE
#
# PSI4: an ab initio quantum chemistry software package
#
# This program is free software; you can redistribute it and/or modify
# it under the terms of the GNU General Public License as published by
# the Free Software Foundation; either version 2 of the License, or
# (at your option) any later version.
#
# This program is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
# GNU General Public License for more details.
#
# You should have received a copy of the GNU General Public License along
# with this program; if not, write to the Free Software Foundation, Inc.,
# 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
#
#@END LICENSE
#

import psi4
import p4const
from math import sqrt, pi
from diatomic_fits import *

[docs]def anharmonicity(rvals, energies, mol = None): """Generates spectroscopic constants for a diatomic molecules. Fits a diatomic potential energy curve using either a 5 or 9 point Legendre fit, locates the minimum energy point, and then applies second order vibrational perturbation theory to obtain spectroscopic constants. The r values provided must bracket the minimum energy point, or an error will result. A dictionary with the following keys, which correspond to spectroscopic constants, is returned: :type rvals: list :param rvals: The bond lengths (in Angstrom) for which energies are provided of length either 5 or 9 but must be the same length as the energies array :type energies: list :param energies: The energies (Eh) computed at the bond lengths in the rvals list :returns: (*dict*) Keys: "re", "r0", "we", "wexe", "nu", "ZPVE(harmonic)", "ZPVE(anharmonic)", "Be", "B0", "ae", "De" corresponding to the spectroscopic constants in cm-1 """ angstrom_to_bohr = 1.0 / p4const.psi_bohr2angstroms angstrom_to_meter = 10e-10; if len(rvals) != len(energies): raise Exception("The number of energies must match the number of distances") npoints = len(rvals) if npoints != 5 and npoints != 9: raise Exception("Only 5- or 9-point fits are implemented right now") psi4.print_out("\n\nPerforming a %d-point fit\n" % npoints) psi4.print_out("\nOptimizing geometry based on current surface:\n\n"); if (npoints == 5): optx = rvals[2] elif (npoints == 9): optx = rvals[4] # Molecule can be passed in be user. Look at the function definition above. if mol == None: mol = psi4.get_active_molecule() natoms = mol.natom() if natoms != 2: raise Exception("The current molecule must be a diatomic for this code to work!") m1 = mol.mass(0) m2 = mol.mass(1) maxit = 30 thres = 1.0e-9 for i in range(maxit): if (npoints == 5): grad= first_deriv_5pt(rvals, energies, optx) secd = second_deriv_5pt(rvals, energies, optx) energy = function_5pt(rvals, energies, optx) elif (npoints == 9): grad = first_deriv_9pt(rvals, energies, optx) secd = second_deriv_9pt(rvals, energies, optx) energy = function_9pt(rvals, energies, optx) psi4.print_out(" E = %20.14f, x = %14.7f, grad = %20.14f\n" % (energy, optx, grad)) if abs(grad) < thres: break optx -= grad / secd; psi4.print_out(" Final E = %20.14f, x = %14.7f, grad = %20.14f\n" % (function_5pt(rvals, energies, optx), optx, grad)); if optx < min(rvals): raise Exception("Minimum energy point is outside range of points provided. Use a lower range of r values.") if optx > max(rvals): raise Exception("Minimum energy point is outside range of points provided. Use a higher range of r values.") if (npoints == 5): energy = function_5pt(rvals, energies, optx) first = first_deriv_5pt(rvals, energies, optx) secd = second_deriv_5pt(rvals, energies, optx) * p4const.psi_hartree2aJ third = third_deriv_5pt(rvals, energies, optx) * p4const.psi_hartree2aJ fourth = fourth_deriv_5pt(rvals, energies, optx) * p4const.psi_hartree2aJ elif (npoints == 9): energy = function_9pt(rvals, energies, optx) first = first_deriv_9pt(rvals, energies, optx) secd = second_deriv_9pt(rvals, energies, optx) * p4const.psi_hartree2aJ third = third_deriv_9pt(rvals, energies, optx) * p4const.psi_hartree2aJ fourth = fourth_deriv_9pt(rvals, energies, optx) * p4const.psi_hartree2aJ psi4.print_out("\nEquilibrium Energy %20.14f Hartrees\n" % energy) psi4.print_out("Gradient %20.14f\n" % first) psi4.print_out("Quadratic Force Constant %14.7f MDYNE/A\n" % secd) psi4.print_out("Cubic Force Constant %14.7f MDYNE/A**2\n" % third) psi4.print_out("Quartic Force Constant %14.7f MDYNE/A**3\n" % fourth) hbar = p4const.psi_h / (2.0 * pi) mu = ((m1*m2)/(m1+m2))*p4const.psi_amu2kg we = 5.3088375e-11*sqrt(secd/mu) wexe = (1.2415491e-6)*(we/secd)**2 * ((5.0*third*third)/(3.0*secd)-fourth) # Rotational constant: Be I = ((m1*m2)/(m1+m2)) * p4const.psi_amu2kg * (optx * angstrom_to_meter)**2 B = p4const.psi_h / (8.0 * pi**2 * p4const.psi_c * I) # alpha_e and quartic centrifugal distortion constant ae = -(6.0 * B**2 / we) * ((1.05052209e-3*we*third)/(sqrt(B * secd**3))+1.0) de = 4.0*B**3 / we**2 # B0 and r0 (plus re check using Be) B0 = B - ae / 2.0 r0 = sqrt(p4const.psi_h / (8.0 * pi**2 * mu * p4const.psi_c * B0)) recheck = sqrt(p4const.psi_h / (8.0 * pi**2 * mu * p4const.psi_c * B)) r0 /= angstrom_to_meter; recheck /= angstrom_to_meter; # Fundamental frequency nu nu = we - 2.0 * wexe; zpve_nu = 0.5 * we - 0.25 * wexe; psi4.print_out("\nre = %10.6f A check: %10.6f\n" % (optx, recheck)) psi4.print_out("r0 = %10.6f A\n" % r0) psi4.print_out("we = %10.4f cm-1\n" % we) psi4.print_out("wexe = %10.4f cm-1\n" % wexe) psi4.print_out("nu = %10.4f cm-1\n" % nu) psi4.print_out("ZPVE(nu) = %10.4f cm-1\n" % zpve_nu) psi4.print_out("Be = %10.4f cm-1\n" % B) psi4.print_out("B0 = %10.4f cm-1\n" % B0) psi4.print_out("ae = %10.4f cm-1\n" % ae) psi4.print_out("De = %10.7f cm-1\n" % de) results = { "re" : optx, "r0" : r0, "we" : we, "wexe" : wexe, "nu" : nu, "ZPVE(harmonic)" : zpve_nu, "ZPVE(anharmonic)" : zpve_nu, "Be" : B, "B0" : B0, "ae" : ae, "De" : de } return results