Question about force-gradient

[abaillod]

Hello,

I am implementing the analytical derivatives for evaluating the force-gradient when using Sophia'a representation. To test it, I am comparing the analytical derivatives to finite differences, using Lcheck=6.

I remember doing the same when implementing the force gradient for the current constraint (Lconstraint=3). I had to enforce the magnetic axis to stay fixed while the finite differences were evaluated - otherwise there would be a miss-match between the analytical and the finit-difference derivatives.

I realized that this is not true in the master branch anymore - to get an agreement, the coordinate axis has to be allowed to move. This change occurred between v3.00 and v3.10 - see the following lines:

https://github.com/PrincetonUniversity/SPEC/blob/66665e217bf602da8e6241694e6b6065cf0f86cc/src/dforce.f90#L781-L786

I have been able to track the change to the commit ce0a8040fd83e59354410c0ea24829adce1c89dc, made by Zhisong.

My question is: why? What is the "correct" force gradient now; the finite differences with moving coordinate axis?

[zhisong]

This is because the axis is determined completely by the first interface. If you take the derivative w.r.t. the first interface, that applies to the axis as well.

So yes, you need to take finite difference with moving coordinate axis.

[abaillod]

Sure; I don't get however why, when implementing the force gradient for the current constraint I had to keep the axis fixed, see these lines, from the release 3.00:

https://github.com/PrincetonUniversity/SPEC/blob/e63128ca34998ba2f683538ae166fddf9e480ba5/dforce.f90#L708-L712

I guess there has been a change in the logic of SPEC, and the coordinate axis is now re-evaluated at every force calculation while it was not necessarily the case when I was working on the current constraint.

[zhisong]

You are correct. This feature was newly added.

[abaillod]

Thanks! I close this thread then.