Add default for properties origin at the COC

[jaclark5]

Description

Make SCF properties calculated from the COC (center of charge) of the molecule rather than the origin. This isn't an issue for dipole moments of neutral molecules, but unexpectedly results in very large SCF multipole moments.

I understand if this should be an "us" problem and we should remember to add the keyword, but this is very unexpected behevior. Both Gaussian and ORCA report multiple moments from the center of nuclear charge, so I think this should be consistent.

Changelog description

Set psi4 properties_origin to "center of charge" instead of origin.

Status

• Code base linted
• Ready to go

[amcisaac]

This has bitten me a number of times too, especially in conjunction with QCElemental's default of fix_com False for molecules and inputs, which (if I understand correctly) is not the default for Psi4 and makes the behavior especially unexpected. I'd support a change to the default, although of course it's a very wide-reaching change. If nothing else it would be nice to get a warning or something

[jaclark5]

This change has been verified to work in DS: 464, the following records have completed: [147472201, 147472680, 147472714, 147472946, 147472929, 147472957, 147473071, 147472975, 147472963, 147473206, 147473178, 147473494, 147473396, 147473319, 147473394, 147473291]

[loriab]

Talking with Ben, I think I understand the problem, and yes properties at the origin is weird with fix_com/fix_orientation. (Psi4's default is to reorient to a standard frame, so the origin makes more sense there.) Is center-of-charge much preferred over center-of-mass? I know COC would have been a better choice to start with for isotopic substitution, but everything else in Psi4 is COM-based.

[jaclark5]

@loriab I agree that COM makes more sense, I suggest COC because it's what the other packages do so psi4 would be consistent. We discussed in an internal meeting and Bill Swope suggested that COC is chosen because it's independent of isotope.

[codecov]

Codecov Report

❌ Patch coverage is 26.47059% with 25 lines in your changes missing coverage. Please review.
✅ Project coverage is 30.42%. Comparing base (0a489b2) to head (cf9c666).
⚠️ Report is 1 commits behind head on master.

🚀 New features to boost your workflow:

• ❄️ Test Analytics: Detect flaky tests, report on failures, and find test suite problems.

[loriab]

Hi, I've rebased this, standardized it on QCSchema v1, and started another test.

I definitely agree with the problem -- you provided a geometry, didn't ask for it to have its origin respected, and yet the molecule is frozen at that geometry, and properties are calculated respecting the input origin. The difficulty is that the "you didn't ask for it to have its origin respected" is imposed by QCFractal (here, I think https://github.com/MolSSI/QCFractal/blob/main/qcfractal/qcfractal/components/molecules/db_models.py#L56), not QCElemental or QCEngine. Below shows an ordinary molecule input that gets adjusted to COM and thus the properties calc from (0,0,0) makes sense. From QCEngine's POV, the PR is saying even if the user specifically asks that the origin be respected (fix_*=T), ignore that, and compute properties from NUCLEAR_CHARGE (or COM), and that doesn't seem very intuitive. I think this needs solving, but it's probably bigger than just here. If you think I've misanalyzed something, please let me know.

>>> import qcelemental as qcel
>>> import qcengine as qcng
>>> aa = qcel.models.v2.Molecule.from_data("He 0, 0, 0\nNe 1, 0, 0")
>>> bb = qcel.models.v2.AtomicInput(**{
... "molecule": aa,
... "specification": {
...     "driver": "energy",
...     "model": {"method": "hf", "basis": "cc-pvdz"},
... },
... })
>>> atres = qcng.compute(bb, "psi4")
>>> print(atres.stdout)
...
  ==> Geometry <==

    Molecular point group: c2v
    Full point group: C_inf_v

    Geometry (in Bohr), charge = 0, multiplicity = 1:

       Center              X                  Y                   Z               Mass       
    ------------   -----------------  -----------------  -----------------  -----------------
         HE           0.000000000000     0.000000000000    -1.574501697116     4.002603254130
         NE           0.000000000000     0.000000000000     0.315224432884    19.992440176200

  Running in c2v symmetry.

  Rotational constants: A = ************  B =      5.05487  C =      5.05487 [cm^-1]
  Rotational constants: A = ************  B = 151541.08358  C = 151541.08358 [MHz]
  Nuclear repulsion =   10.583544187961248

...

  @DF-RHF Final Energy:  -130.93419582222612

   => Energetics <=

    Nuclear Repulsion Energy =             10.5835441879612482
    One-Electron Energy =                -206.6264580667780422
    Two-Electron Energy =                  65.1087180565906891
    Total Energy =                       -130.9341958222261155

Computation Completed


Properties will be evaluated at   0.000000,   0.000000,   0.000000 [a0]

Properties computed using the SCF density matrix


 Multipole Moments:

 ------------------------------------------------------------------------------------
     Multipole            Electronic (a.u.)      Nuclear  (a.u.)        Total (a.u.)
 ------------------------------------------------------------------------------------

 L = 1.  Multiply by 2.5417464519 to convert [e a0] to [Debye]
 Dipole X            :          0.0000000            0.0000000            0.0000000
 Dipole Y            :          0.0000000            0.0000000            0.0000000
 Dipole Z            :          0.2402701            0.0032409            0.2435110
 Magnitude           :                                                    0.2435110

 ------------------------------------------------------------------------------------

[jaclark5]

@loriab let me lay out what I understand and get @bennybp opinion.

We agree that it's a problem that the origin is used to calculate the scf properties by default. This could be resolved / changed in a few places:

1. Downstream: QCFractal: It sounds like you've chased down the most down stream source of why calculations on QCArchive calculate properties from the origin. It's because in the QCFractal db_model.MoleculeORM the class variable is set fix_com = column_property(True). One might then think that this default could be changed, but this setting applied to all molecules in the entire database, @bennybp would that be feasible? It does seem short weird that QCFractal sets fix_com=True no matter what for all molecules when this QCElemental property defaults to False ordinarily. This keyword passes through QCEngine and straight to psi4.

2. Upstream: psi4: I've also made a PR about this in Pis4. In Psi4 it is clear that the default center for property calculations is the coordinate origin but can be adjusted to be the COM or COC with the keywords PROPERTIES_ORIGIN["NUCLEAR_CHARGE"]. In other conversations with @bennybp (if I understood correctly) he said that this should not be a psi4 issue since the structure should automatically be translated to the COM as the origin before the calculation. However, this appears to be dependent on the fix_com setting which QCFractal is restraining to True (which prevents translating the structure for the scf center at COM). So we have two ways of controlling it the center for SCF property calculations, PROPERTIES_ORIGIN["NUCLEAR_CHARGE"] and fix_com.

3. Midstream: QCEngine: I anticipated that a change in defaults in psi4 from PROPERTIES_ORIGIN[(0, 0, 0)] to PROPERTIES_ORIGIN["NUCLEAR_CHARGE"] would be an up-hill battle, and at the time thought it was the only way to correct this problem (thanks @loriab for pointing out fix_com). In my tests for this PR I have confirmed that my proposed QCEngine default setting with PROPERTIES_ORIGIN could be overwritten by specifying it to something else and so isn't a permanent setting. However the fact that my tests with PROPERTIES_ORIGIN have been successful shows that it will always override fix_com making that keyword ineffective. By moving forward with this PR, QCEngine would be entering a conflict that it's ignorant of at the moment.

I have confirmed that setting PROPERTIES_ORIGIN["NUCLEAR_CHARGE"] or PROPERTIES_ORIGIN["COM"] in a QCPortal entry will have the desired effect (mentioned in the comments above). I've previously proposed changing the default PROPERTIES_ORIGIN in psi4 or setting a different default in QCEngine, but the former would be disruptive to the community (although align with similar packages) and the later would render the fix_com ineffective. It seems the most straightforward choice would be to change the QCFractal default fix_com=False to align with user expectations, and/or allow it to be changed. If I can get the all clear from @bennybp I'll make a QCFractal PR.

[loriab]

@jaclark5, I agree with your analysis above. iirc, the motivation for fix_=T imposed by qcfractal was to provide a little more standardization among calcs living in the qca database. For example, gradients should be the same btwn qc programs if run as fix_=T; whereas, each program has its own std orientation so unlikely to match at fix_*=F.

I think the longer-term sol'n is:
(1) change psi4's std orientation overall to use center of charge, rather than center of mass, since it's a cleaner approach for isotope purposes. I brought this up at PsiCon last month and no one had strong opinions. Another complication to work through that I realized ysterday is that for the purposes of std orientation and properties_origin, psi4 takes ghosts to have mass. this means that the molecule stays put when you change ghostiness for BSSE counterpoise corrections. But for properties_origin, ghosts make a difference for COM (to match std orient) but don't make a difference for NUCLEAR_CHARGES. Need to think more on if that's the right decision.
(2) bring properties into the fold of what's standardized by fix_* commands in QCEngine. That gradient assertion I made above gets checked with test_align, but properties don't right now.

[jaclark5]

@loriab Oh I'm glad to hear that you floated the idea of COC at PsiCon, I hope that went well!

I'm not sure what you mean by the second point, but it sounds like an issue should be made in this repo for it, and that this current PR should be closed and not pursued.

[jaclark5]

Will not pursue issue in this repo