An atomic-orbital based Lagrangian approach for calculating geometric gradients of linear response properties