Closed cortner closed 3 months ago
Comments carried over from the AC PR:
(sys, calc)
instead of (calc, sys)
(to discuss) Main comment from @tjjarvinen to be discussed:
For main point the PR is ok. I marked some comments to the code.
But there is a big issue with units.
Basicly you now implicitely expect that force unit is energy unit/position unit of AtomsBase structure. This leads to a issue where the results change, if you just change the position unit for AtomsBase system.
In my opinnion units should not be stripped. It just leads to errors and because AtomsCalculators and AtomsBase have been designed units in mind, it also just complicates the code.
Here is a example of this effect
using ORCAcalculator
using AtomsBase
using AtomsCalculators
using Unitful
using UnitfulAtomic
using AtomsCalculators.AtomsCalculatorsTesting
ox = ORCAexecutable()
om = ORCAmethod("blyp def2-TZVPP TIGHTSCF defgrid3")
h1 = isolated_system([
:H => [0, 0, 0.]u"Å",
:H => [0, 0, 1.]u"Å"
])
h2 = isolated_system([
:H => auconvert.([0, 0, 0.]u"Å"),
:H => auconvert.([0, 0, 1.]u"Å")
])
orca = ORCAcalculator.ORCAcalculatorbundle(ox, om)
# These two are different
AtomsCalculatorsTesting.fdtest(orca, h1; test_virial=false)
AtomsCalculatorsTesting.fdtest(orca, h2; test_virial=false)
In my opinnion units should not be stripped
true, this needs to be fixed if they actually are stripped? But where does this actually happen? The core FD test cannot handle units and shouldn't have to handle units. But the wrappers should manage the units, which they do. So I'm not yet clear where exactly the problem is.
Basicly you now implicitely expect that force unit is energy unit/position unit of AtomsBase structure
I think this is the correct procedure. The calculator gets to pick the energy unit, but the system has already decided on the force unit. Nothing to do about it.
In my view the calculator is responsible to either return the forces int he correct unit or to throw an error.
This is part of the bigger discussion how to sort out units. Happy to hear your thoughts.
We can implement a calculator to wrap other calculators and conversts their output to predefined units. This would be fast to implement and solves your issue here.
I don't understand the issue yet. Let's hold off on this please?
Can I just check to be sure; your concern is that the calculator might return a different unit from uE / uL where uL is the system length unit? And that somehow leads to problems with the FD test. (which isn't yet clear to me, but I can think it through again)
I believe the last two commits fix the unit conversion issue - at least in principle. Obviously there can still be bugs. Please review and merge if you are happy.
In the meantime I will merge this into the sitepotentail branch so I continue there.
There is a compatibility test. It has to be somewhere. The system and calculator must have the same units. We can discuss where that test should be. Maybe inside the calculator
My opinion is that we should not start to split AtomsCalculators features, so that they only work with some calculators. Thus, I am against merging anything like that into a general package line AtomsCalculatorsUtilities.
edit. rephrased comment
I violently disagree with that statement. Please think about what you are saying.
Most of the utilities are not for general calculators. And the general philosophy of Julia is that things work whenever the suitable conditions are met/ eg methods implemented.
To calculate FD gradient you only need potential energy calculation implemented. Every calculator does that. Thus I see no reason why we could not get it to work with a general calculator.
edit. rephrased comment.
Regarding your ORCA example: the @assert that fails is actually just an assert for documentation. The units are converted to match:
uL_sys = unit(position(sys, 1)[1])
h0 = uconvert(uL_sys, h0)
~I suspect~ This is another infinity bug and it has nothing to do with units.
These infinities are an abomination. This is now fixed, but yet another ad hoc hack. This needs to be removed from the AtomsBase interface.
With the current version using
using ORCAcalculator
using AtomsBase
using AtomsCalculators
using Unitful
using UnitfulAtomic
using AtomsCalculatorsUtilities.Testing
ox = ORCAexecutable()
om = ORCAmethod("blyp def2-TZVPP TIGHTSCF defgrid3")
h1 = isolated_system([
:H => [0, 0, 0.]u"Å",
:H => [0, 0, 1.]u"Å"
])
h2 = isolated_system([
:H => auconvert.([0, 0, 0.]u"Å"),
:H => auconvert.([0, 0, 1.]u"Å")
])
orca = ORCAcalculator.ORCAcalculatorbundle(ox, om)
Testing.fdtest(h1, orca; test_virial=false)
I get
Forces finite-difference test
---------|-----------
h | error
---------|-----------
[file orca_tools/Tool-Misc/qccc2ic.cpp, line 654]: Zero distance between atoms 2 and 1 in Cartesian2Internal
ERROR: failed process: Process(`/home/teemu/.local/opt/orca-5.0.4/orca /tmp/jl_V4xD0Q/orca_calculation.inp`, ProcessExited(18)) [18]
Is this related to the infinity bug?
I don't know what causes this. If you want me to debug it then I need a MWE.
Btw, the FD tests need not be universal at this point. This example is not enough to prevent merging. They are useful now as they are.
MWE is basicly the code sniplet I posted. You need to have ORCAcalculator installed ofcourse.
It is not a MWE if I need to install ORCA first and then figure out how to use your ORCA calculator.
That is true. I will try to find out some other way to get to that error.
I have a MWE to reproduce this and will try to fix it.
done - try your test again.
All works now!
The results are now consisten regardless of system unit type and no other errors. That documentation build issue can fixed later (need to add entry to @autodocs
in docs/make.jl
)
I am more than happy to merge this.
thank you. I will also write some docs. But I'd be grateful if you can get started on that. I.e. give them a framework and working docs build. Then I can add a bit.
This PR transfers the FD tests for AtomsCalculators to ACU.