from_zmat: method
from_pybel: method
from_file: method

__init__(self, atoms, coords, bonds=None, obmol=None, charge=None, name=None, zmatrix=None, dipole_surface=None, potential_surface=None, potential_derivatives=None, source_file=None, guess_bonds=True, **kw): 

__repr__(self): 

@property
num_atoms(self): 

@property
atoms(self): 

@property
masses(self): 

@property
bonds(self): 

@property
coords(self): 

@property
sys(self): 

@property
formula(self): 

@property
multiconfig(self): 

@property
name(self): 

@property
force_constants(self): 

@property
potential_derivatives(self): 

@property
potential_surface(self): 

@property
dipole_surface(self): 

@property
center_of_mass(self): 

@property
inertial_axes(self): 

@property
internal_coordinates(self): 

@property
normal_modes(self): 

take_submolecule(self, spec): 

@property
shape(self): 

__len__(self): 

__iter__(self): 

__getitem__(self, item): 

copy(self): 

prop(self, name, *args, **kwargs): 

load_force_constants(self, file=None): 

load_potential_derivatives(self, file=None): 

load_normal_modes(self, file=None): 

load_potential_surface(self): 

load_dipole_surface(self): 

principle_axis_frame(self, sel=None, inverse=False): 

eckart_frame(self, mol, sel=None, inverse=False): 

get_embedded_molecule(self, ref=None): 

plot(self, *geometries, figure=None, bond_radius=0.1, atom_radius_scaling=0.25, atom_style=None, bond_style=None, mode='fast', objects=False, **plot_ops): 

get_normal_modes(self, **kwargs): 

test_log_water = TestManager.test_data("water_OH_scan.log")
test_log_freq = TestManager.test_data("water_freq.log")
test_HOD = TestManager.test_data("HOD_freq.fchk")
test_fchk = TestManager.test_data("water_freq.fchk")
test_log_h2 = TestManager.test_data("outer_H2_scan_new.log")

def test_ImportMolecule(self):

    n = 3 # water
    m = Molecule.from_file(self.test_fchk)
    self.assertEquals(m.atoms, ("O", "H", "H"))

def test_Eckart(self):
    scan_file = TestManager.test_data("tbhp_030.log")
    ref_file = TestManager.test_data("tbhp_180.fchk")

    scan = Molecule.from_file(scan_file)
    ref = Molecule.from_file(ref_file)
    sel = np.where(ref.masses > 3)[0]

    pax_rot = ref.principle_axis_frame(sel=sel) #type: MolecularTransformation
    rot_ref = pax_rot.apply(ref)

    #transf = scan.principle_axis_frame(sel=sel)
    transf = scan.eckart_frame(ref, sel=sel)
    tf_test = transf[0].transformation_function

    tf_mat = tf_test.transform
    self.assertTrue(np.allclose(tf_mat@tf_mat.T - np.eye(3), 0.))
    self.assertEquals(tf_test.transf.shape, (4, 4))

    for t, m in zip(transf, scan):
        new_mol = t(m)
        # rot_ref.guess_bonds = False
        # ref.guess_bonds = False
        # m.guess_bonds = False
        # new_mol.guess_bonds = False
        # m = m #type: Molecule
        # g1, a, b = ref.plot()
        # # ref.plot(figure=g1)
        # rot_ref.plot(figure=g1)
        # g, a, b = new_mol.plot()
        # rot_ref.plot(figure=g, atom_style=dict(color='black'))
        # g.show()
        fuckup = np.linalg.norm(new_mol.coords[sel] - rot_ref.coords[sel])
        self.assertLess(fuckup/len(sel), .1)

def test_EckartEmbedDipoles(self):
    scan_file = TestManager.test_data("tbhp_030.log")
    ref_file = TestManager.test_data("tbhp_180.fchk")

    scan = Molecule.from_file(scan_file)
    ref = Molecule.from_file(ref_file)
    sel = np.where(ref.masses>3)[0]
    pax_rot = ref.principle_axis_frame(sel=sel, inverse=True)  # type: MolecularTransformation
    rot_ref = pax_rot.apply(ref)

    transf = scan.eckart_frame(rot_ref, sel=sel)

    carts, dips = DipoleSurface.get_log_values(scan_file, keys=("StandardCartesianCoordinates", "OptimizedDipoleMoments"))
    rot_dips = np.array([ np.dot(t.transformation_function.transform, d) for t,d in zip(transf, dips) ])
    self.assertTrue(np.allclose(np.linalg.norm(dips, axis=1)-np.linalg.norm(rot_dips, axis=1), 0.))

    ### Visualize dipole surface
    dists = np.linalg.norm(carts[1:, 5] - carts[1:, 6], axis=1)
    Graphics.default_style['image_size'] = 575
    g = GraphicsGrid(nrows=1, ncols=2, padding=((.075, 0), (0, .45)))
    p = Plot(dists, rot_dips[:, 0], figure=g[0, 0])
    Plot(dists, rot_dips[:, 1], figure=p)
    Plot(dists, rot_dips[:, 2], figure=p)
    p2= Plot(dists, dips[:, 0], figure=g[0, 1])
    Plot(dists, dips[:, 1], figure=p2)
    Plot(dists, dips[:, 2], figure=p2)
    g.show()

def test_NormalModes(self):
    n = 3 # water
    with GaussianFChkReader(self.test_fchk) as reader:
        parse = reader.parse(("ForceConstants", "AtomicMasses"))
    fcs = parse["ForceConstants"].array
    masses = parse["AtomicMasses"]

    nms = MolecularNormalModes.from_force_constants(fcs, masses)
    # ArrayPlot(nms.matrix).show()
    self.assertEquals(nms.matrix.shape, (3*n, 3*n))
    
def test_Plotting(self):

    g = Graphics3D(
        image_size=[1500, 1500],
        plot_range=[[-10, 10]]*3,
        backend="VTK"
        )
    h5 = Molecule.from_file(
        self.test_log_h2,
        # self.test_fchk,
        # bonds = [
        #     [0, 1, 1],
        #     [0, 2, 1]
        # ]
    )
    h5.plot(
        figure=g
        # mode='3D',
        # bond_style= { "circle_points": 24 },
        # atom_style= { "sphere_points": 24 }
    )
    m = Molecule.from_file(
        self.test_fchk,
        bonds = [
            [0, 1, 1],
            [0, 2, 1]
        ]
    )
    m.plot(
        figure=g
        # mode='3D',
        # bond_style= { "circle_points": 24 },
        # atom_style= { "sphere_points": 24 }
        )
    # g.show()

def test_BondGuessing(self):
    m = Molecule.from_file(self.test_fchk)
    self.assertEquals(m.bonds, [[0, 1, 1], [0, 2, 1]])

def test_Frags(self):
    m = Molecule.from_file(self.test_fchk)
    self.assertEquals(len(m.prop("fragments")), 1)

def test_HODModes(self):
    # oops fucked up getting D out
    m = Molecule.from_file(self.test_HOD, bonds=[[0, 1, 1], [0, 2, 1]])
    modes = m.normal_modes
    self.assertEquals(m.atoms, ("O", "H", "D"))
    self.assertEquals(tuple(np.round(modes.freqs)), (1422.0, 2810.0, 3874.0))

def test_H2OModes(self):
    m = Molecule.from_file(self.test_fchk, bonds=[[0, 1, 1], [0, 2, 1]])
    modes = m.normal_modes
    self.assertEquals(m.atoms, ("O", "H", "H"))
    self.assertEquals(tuple(np.round(modes.freqs)), (1622.0, 3803.0, 3938.0))

def test_RenormalizeGaussianModes(self):

    with GaussianFChkReader(self.test_HOD) as gr:
        parse = gr.parse(["Coordinates", "Gradient", "AtomicMasses",
                          "ForceConstants", "ForceDerivatives", "VibrationalModes", "VibrationalData"])

    coords = UnitsData.convert("Angstroms", "AtomicUnitOfLength") * parse["Coordinates"]
    masses = UnitsData.convert("AtomicMassUnits", "AtomicUnitOfMass") * parse["AtomicMasses"]
    modes = parse["VibrationalModes"].T
    freqs = parse["VibrationalData"]["Frequencies"]
    fcs = parse["ForceConstants"].array
    sad = UnitsData.convert("Hartrees", "Wavenumbers") * np.sqrt(np.diag(np.dot(np.dot(modes.T, fcs), modes)))
    modes = modes * freqs/sad
    print( UnitsData.convert("Hartrees", "Wavenumbers") * np.sqrt(np.diag(np.dot(np.dot(modes.T, fcs), modes))))

    masses = np.broadcast_to(masses, (len(masses), 3)).T.flatten()
    # print(modes-np.linalg.pinv(modes).T)
    print(np.dot(np.dot(modes.T, np.diag(masses)), modes))

    modes_2 = Molecule.from_file(self.test_HOD).get_normal_modes(normalize=False)
    mm = modes_2._basis.matrix

    print(np.dot(np.dot(mm.T, np.diag(masses)), mm))
    print(UnitsData.convert("Hartrees", "Wavenumbers") * np.sqrt(np.diag(np.dot(np.dot(mm.T, fcs), mm))))
    # print(modes._basis.matrix.T.dot(m.force_constants).shape)

def test_VisualizeNormalModes(self):

    from Psience.Molecools.Vibrations import MolecularVibrations, MolecularNormalModes
    from McUtils.Plots import GraphicsGrid, Graphics3D

    m = Molecule.from_file(self.test_fchk, bonds = [[0, 1, 1], [0, 2, 1]])

    with GaussianFChkReader(self.test_fchk) as reader:
        parse = reader.parse(("VibrationalModes", "VibrationalData"))
    modes = parse["VibrationalModes"].T

    test_freqs = parse["VibrationalData"]["Frequencies"]

    nms = m.normal_modes
    realvibs = MolecularVibrations(m, basis=MolecularNormalModes(modes, freqs=test_freqs))

    plot_vibrations = False
    if plot_vibrations:
        nmodes = 1
        mode_start = 0
        g = GraphicsGrid(nrows=2, ncols=nmodes,
                         graphics_class=Graphics3D,
                         plot_range = [[-2, 2], [-2, 2], [-2, 2]],
                         fig_kw = dict(figsize = (17, 5)),
                         tighten = True
                         )

        for i in range(nmodes):
            nms.visualize(step_size=.1, figure = g[0, i], which=mode_start + i,
                          anim_opts= dict(interval = 10)
                          )

        for i in range(nmodes):
            realvibs.visualize(step_size=.1, figure = g[1, i], which= mode_start+i,
                               anim_opts= dict(interval = 10)
                               )

        g.show()

    self.assertEquals(tuple(round(a, 4) for a in nms.freqs), tuple(round(a, 4) for a in test_freqs))