#include "lbfgs.h"
#include "source_pw/module_pwdft/global.h"
#include "source_base/matrix3.h"
#include "source_io/module_parameter/parameter.h"
#include "ions_move_basic.h"
#include "source_cell/update_cell.h"
#include "source_cell/print_cell.h" // mohan add 2025-06-19  

void LBFGS::allocate(const int _size) // initialize H0Hpos0force0force
{
    alpha=70;//default value in ase is 70
    maxstep=PARAM.inp.relax_bfgs_rmax;
    size=_size;
    memory=100;
    iteration=0;
    H = std::vector(3*size, std::vector(3*size, 0.0));
    H0=1/alpha;
    pos = std::vector (size, std::vector(3, 0.0)); 
    pos0 = std::vector(3*size, 0.0);
    pos_taud = std::vector (size, std::vector(3, 0.0)); 
    pos_taud0 = std::vector(3*size, 0.0);
    dpos = std::vector(size, std::vector(3, 0.0));
    force0 = std::vector(3*size, 0.0);
    force = std::vector(size, std::vector(3, 0.0));
    steplength = std::vector(size, 0.0);
    //l_search.init_line_search();
}

void LBFGS::relax_step(const ModuleBase::matrix _force,UnitCell& ucell,const double &etot)

{
    get_pos(ucell,pos);  
    get_pos_taud(ucell,pos_taud);
    //solver=p_esolver;
    ucell.ionic_position_updated = true;
    for(int i = 0; i < _force.nr; i++)
    {
        for(int j=0;jupdate_pos(ucell);
    this->calculate_largest_grad(_force,ucell);
    this->is_restrain(dpos);  
    // mohan add 2025-06-22
    unitcell::print_tau(ucell.atoms,ucell.Coordinate,ucell.ntype,ucell.lat0,GlobalV::ofs_running);
}

void LBFGS::get_pos(UnitCell& ucell,std::vector& pos)
{
    int k=0;
    for(int i=0;i= maxstep)
    {
        double scale = maxstep / a;
        for(int i = 0; i < size; i++)
        {
            for(int j=0;j