Ensure that molecule orientation remains fixed for MDI

taylor-a-barnes commented 10 months ago

Description

This is a simple bug fix to ensure that molecules remain fixed during MDI calculations. Otherwise, there is a possibility that certain codes may rotate the molecule, which renders quantities such as forces unusable by external drivers.

Changelog description

Ensure that molecule orientation remains fixed for MDI.

Status

[x] Code base linted
[x] Ready to go

codecov[bot] commented 10 months ago

Codecov Report

Merging #428 (f8ae2a0) into master (cbec059) will decrease coverage by 1.58%. Report is 1 commits behind head on master. The diff coverage is 80.00%.

Additional details and impacted files

loriab commented 10 months ago

Hi @taylor-a-barnes, you will have noticed that this change causes the mdi/psi4/qcng/gradient test to fail, and I can reproduce that locally. After the tests here https://github.com/MolSSI/QCEngine/blob/master/qcengine/programs/tests/test_alignment.py#L57-L60 I'm fairly confident that psi4 and qcng are working as expected for fix_* fields. Could psi4's orientation handling not match what you're expecting? From the psi4/mdi (no qcng) test case https://github.com/psi4/psi4/blob/master/tests/pytests/test_mdi.py#L13 I put together a little demo of how the forces are returned from some input waters of different setup. Anything jarring?

import psi4

psi4.set_output_file("forces_output.dat")
psi4.set_options({
    'reference': 'uhf',
    'scf_type': 'direct',
    'e_convergence': 1e-12,
})
scf_method = "scf/cc-pVDZ"

#####

smol = """
    0 1
    O
    H 1 1.0
    H 1 1.0 2 104.5
"""
mol_zmat = psi4.geometry(smol)

g = psi4.gradient(scf_method, molecule=mol_zmat)
print("\n<<<  ZMAT (from initial mol)  >>>")
print(smol)
print(f"{mol_zmat.nuclear_repulsion_energy()=}")
print(psi4.nppp(mol_zmat.geometry().np))
print(psi4.nppp(g.np))

#####

smol = """
    units au
    O  0.0,  0.0,                -0.12947688966627896
    H  0.0, -1.4941867446308548,  1.027446098871551
    H  0.0,  1.4941867446308548,  1.027446098871551
"""
mol_cart_4z = psi4.geometry(smol)

g = psi4.gradient(scf_method, molecule=mol_cart_4z)
print("\n<<<  CART 4-zeros (from <COORDS cmd)  >>>")
print(smol)
print(f"{mol_cart_4z.nuclear_repulsion_energy()=}")
print(psi4.nppp(mol_cart_4z.geometry().np))
print(psi4.nppp(g.np))

#####

smol = """
    units au
    O  0.1,  0.0,                -0.12947688966627896
    H  0.0, -1.4341867446308548,  1.027446098871551
    H  0.0,  1.4941867446308548,  1.027446098871551
"""
mol_cart_3z = psi4.geometry(smol)

g = psi4.gradient(scf_method, molecule=mol_cart_3z)
print("\n<<<  CART 3-zeros (from >COORDS cmd)  >>>")
print(smol)
print(f"{mol_cart_3z.nuclear_repulsion_energy()=}")
print(psi4.nppp(mol_cart_3z.geometry().np))
print(psi4.nppp(g.np))

#####

smol = """
    units au
    O  0.0,  0.0,                -0.12947688966627896
    H  0.0, -1.4941867446308548,  1.027446098871551
    H  0.0,  1.4941867446308548,  1.027446098871551
    no_com
    no_reorient
"""
mol_cart_4z_fix = psi4.geometry(smol)

g = psi4.gradient(scf_method, molecule=mol_cart_4z_fix)
print("\n<<<  CART 4-zeros, fixed  >>>")
print(smol)
print(f"{mol_cart_4z_fix.nuclear_repulsion_energy()=}")
print(psi4.nppp(mol_cart_4z_fix.geometry().np))
print(psi4.nppp(g.np))

#####

smol = """
    units au
    O  0.1,  0.0,                -0.12947688966627896
    H  0.0, -1.4341867446308548,  1.027446098871551
    H  0.0,  1.4941867446308548,  1.027446098871551
    no_com
    no_reorient
"""
mol_cart_3z_fix = psi4.geometry(smol)

g = psi4.gradient(scf_method, molecule=mol_cart_3z_fix)
print("\n<<<  CART 3-zeros, fixed  >>>")
print(smol)
print(f"{mol_cart_3z_fix.nuclear_repulsion_energy()=}")
print(psi4.nppp(mol_cart_3z_fix.geometry().np))
print(psi4.nppp(g.np))

<<<  ZMAT (from initial mol)  >>>

    0 1
    O
    H 1 1.0
    H 1 1.0 2 104.5

mol_zmat.nuclear_repulsion_energy()=8.801465564567374
[[ 0.          0.         -0.12947689]
 [ 0.         -1.49418674  1.0274461 ]
 [ 0.          1.49418674  1.0274461 ]]
[[-0.          0.         -0.05958518]
 [ 0.         -0.04304646  0.02979259]
 [-0.          0.04304646  0.02979259]]

<<<  CART 4-zeros (from <COORDS cmd)  >>>

    units au
    O  0.0,  0.0,                -0.12947688966627896
    H  0.0, -1.4941867446308548,  1.027446098871551
    H  0.0,  1.4941867446308548,  1.027446098871551

mol_cart_4z.nuclear_repulsion_energy()=8.801465564567374
[[ 0.          0.         -0.12947689]
 [ 0.         -1.49418674  1.0274461 ]
 [-0.          1.49418674  1.0274461 ]]
[[-0.          0.         -0.05958518]
 [ 0.         -0.04304646  0.02979259]
 [-0.          0.04304646  0.02979259]]

<<<  CART 3-zeros (from >COORDS cmd)  >>>

    units au
    O  0.1,  0.0,                -0.12947688966627896
    H  0.0, -1.4341867446308548,  1.027446098871551
    H  0.0,  1.4941867446308548,  1.027446098871551

mol_cart_3z.nuclear_repulsion_energy()=8.904181021442504
[[ 0.12995966 -0.00335745  0.        ]
 [-1.03127709 -1.43754419  0.        ]
 [-1.03127709  1.4908293   0.        ]]
[[ 0.05195064 -0.01775852  0.        ]
 [-0.01946545 -0.02409327  0.        ]
 [-0.03248519  0.04185179  0.        ]]

<<<  CART 4-zeros, fixed  >>>

    units au
    O  0.0,  0.0,                -0.12947688966627896
    H  0.0, -1.4941867446308548,  1.027446098871551
    H  0.0,  1.4941867446308548,  1.027446098871551
    no_com
    no_reorient

mol_cart_4z_fix.nuclear_repulsion_energy()=8.801465564567374
[[ 0.          0.         -0.12947689]
 [ 0.         -1.49418674  1.0274461 ]
 [-0.          1.49418674  1.0274461 ]]
[[-0.          0.         -0.05958518]
 [ 0.         -0.04304646  0.02979259]
 [-0.          0.04304646  0.02979259]]

<<<  CART 3-zeros, fixed  >>>

    units au
    O  0.1,  0.0,                -0.12947688966627896
    H  0.0, -1.4341867446308548,  1.027446098871551
    H  0.0,  1.4941867446308548,  1.027446098871551
    no_com
    no_reorient

mol_cart_3z_fix.nuclear_repulsion_energy()=8.904181021442504
[[ 0.1         0.         -0.12947689]
 [ 0.         -1.43418674  1.0274461 ]
 [ 0.          1.49418674  1.0274461 ]]
[[ 0.00447373 -0.01775852 -0.05175765]
 [-0.00167627 -0.02409327  0.01939314]
 [-0.00279746  0.04185179  0.03236451]]

taylor-a-barnes commented 10 months ago

@loriab Thanks for looking into this. After running some more tests, I've concluded that the original expected forces in the test were represented in an incorrect reference frame. I've updated the comparison in the test to use the correct values.

MolSSI / QCEngine