Skip to content

Commit

Permalink
...EJB
Browse files Browse the repository at this point in the history
  • Loading branch information
@ebylaska
ebylaska committed Nov 2, 2024
1 parent 466436a commit 94e7146
Show file tree
Hide file tree
Showing 5 changed files with 67 additions and 16 deletions.
40 changes: 37 additions & 3 deletions Nwpw/pspw/lib/molecule/Molecule.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -144,7 +144,7 @@ void Molecule::psi_minimize(double *vall, std::ostream &coutput)
bool continue_inner_loop = true;
for (int l=0; continue_inner_loop && l<=(1+(l2-1)*3); ++l)
{
eorb0 = psi_KS_update(ms,k,2,tole,0.001,vall,orb,&error_out,coutput);
eorb0 = psi_KS_update_orb(ms,k,2,tole,0.001,vall,orb,&error_out,coutput);
if (error_out <= tole)
continue_inner_loop = false; // Exit the inner loop
}
Expand Down Expand Up @@ -235,7 +235,7 @@ void Molecule::psi_get_gradient(double *orb, double *vall, double *Horb)

/********************************************
* *
* Molecule::psi_KS_update *
* Molecule::psi_KS_update_orb *
* *
********************************************/
// This routine performs a KS update on orbital i
Expand Down Expand Up @@ -283,7 +283,7 @@ void Molecule::psi_get_gradient(double *orb, double *vall, double *Horb)
* - `mygrid::cc_pack_daxpy()`
* - `myelectron::get_myke()->ke_precondition()`
*/
double Molecule::psi_KS_update(const int ms, const int k, const int maxit_orb, const double maxerror,
double Molecule::psi_KS_update_orb(const int ms, const int k, const int maxit_orb, const double maxerror,
const double perror, double *vall, double *orb,
double *error_out, std::ostream &coutput)
{
Expand Down Expand Up @@ -403,6 +403,40 @@ double Molecule::psi_KS_update(const int ms, const int k, const int maxit_orb, c
return e0;
}

/********************************************
* *
* Molecule::psi_KS_update *
* *
********************************************/
double Molecule::psi_KS_update(const int maxit_orb, const double maxerror,
const double perror, double *vall,
const int ispin, const int *neq, double *psi,
double *error_out, std::ostream &coutput)

{

for (auto ms=0; ms<ispin; ++ms)
{
int ishift = ms*neq[0]*2*mygrid->PGrid::npack(1);
for (auto i=0; i<neq[ms]; ++i)
{
int indx = 2*mygrid->PGrid::npack(1)*i + ishift;
double *orb = psi + indx;

// orthogonalize to lower orbitals
mygrid->g_project_out_filled_below(psi1, ms, i, orb);

// normalize
double norm = mygrid->cc_pack_dot(1,orb,orb);
norm = 1.0/std::sqrt(norm);
mygrid->c_pack_SMul(1,norm,orb);

double e0 = psi_KS_update_orb(ms, i, maxit_orb, maxerror, perror, vall, orb,
error_out, coutput);
}
}
}

/********************************************
* *
* psi_linesearch_update *
Expand Down
8 changes: 7 additions & 1 deletion Nwpw/pspw/lib/molecule/Molecule.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -121,7 +121,13 @@ class Molecule {

void psi_minimize(double *, std::ostream &);
void psi_get_gradient(double *, double *, double *);
double psi_KS_update(const int, const int, const int, const double, const double, double *, double *, double *, std::ostream &);
double psi_KS_update_orb(const int, const int, const int, const double, const double, double *, double *, double *, std::ostream &);


double psi_KS_update(const int, const double, const double, double *,
const int, const int *, double *,
double *, std::ostream &);


void psi_linesearch_update(double, double, double *, double *, double *, double *);
void psi_sort();
Expand Down
2 changes: 1 addition & 1 deletion Nwpw/pspw/minimizer/cgsd.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -30,7 +30,7 @@ extern double cgsd_bybminimize(Molecule &, Geodesic *, double *, double *,
extern double cgsd_bybminimize2(Molecule &, Geodesic *, double *, double *,
double *, int, int,int,
nwpw_scf_mixing &,
double, double);
double, double, std::ostream &);

} // namespace pwdft
#endif
17 changes: 14 additions & 3 deletions Nwpw/pspw/minimizer/cgsd_bybminimize2.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -31,7 +31,7 @@ double cgsd_bybminimize2(Molecule &mymolecule, Geodesic *mygeodesic, double *E,
double *deltae, double *deltac, int current_iteration,
int ks_it_in, int ks_it_out,
nwpw_scf_mixing &scfmix,
double tole, double tolc)
double tole, double tolc, std::ostream &coutput)
{
bool done = false;
double tmin = 0.0;
Expand All @@ -40,11 +40,19 @@ double cgsd_bybminimize2(Molecule &mymolecule, Geodesic *mygeodesic, double *E,
double sum0, sum1, scale, total_energy;
double dE, max_sigma, min_sigma;
double Eold, dEold, Enew;
double tmin0, deltae0;
double tmin0, deltae0, perror,error_out;

Pneb *mygrid = mymolecule.mygrid;
mygeodesic_ptr = mygeodesic;
double *vall_out;

double ks_deltae = tole;
int ispin = mygrid->ispin;
int *neq = mygrid->neq;
double *vall = mygrid->r_nalloc(ispin);
mymolecule.gen_vall();
mymolecule.get_vall(vall);



// ion-ion energy
double eion = mymolecule.eion();
Expand Down Expand Up @@ -89,6 +97,9 @@ double cgsd_bybminimize2(Molecule &mymolecule, Geodesic *mygeodesic, double *E,
// scf_algorithm,scf_alpha,diis_histories,
// mygrid->ispin,mygrid->n2ft3d,vall_out);

double e0 = mymolecule.psi_KS_update(ks_it_in,ks_deltae,perror,vall,
ispin, neq, mymolecule.psi1,
deltae, coutput);


/* iniitialize blocked cg */
Expand Down
16 changes: 8 additions & 8 deletions Nwpw/pspw/minimizer/cgsd_energy.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -71,7 +71,7 @@ double cgsd_energy(Control2 &control, Molecule &mymolecule, bool doprint, std::o
double kerker_g0 = control.kerker_g0();
int diis_histories = control.diis_histories();
int scf_algorithm = control.scf_algorithm();

double dt = control.time_step();
double dte = dt/sqrt(control.fake_mass());

Expand Down Expand Up @@ -103,8 +103,8 @@ double cgsd_energy(Control2 &control, Molecule &mymolecule, bool doprint, std::o
}

// if (minimizer > 1) pspw_Grsm_list_start()
if ((minimizer == 5) || (minimizer == 8))
it_out = 1;
//if ((minimizer == 5) || (minimizer == 8))
// it_out = 1;

Geodesic12 mygeodesic12(minimizer, &mymolecule, control);

Expand Down Expand Up @@ -378,8 +378,7 @@ double cgsd_energy(Control2 &control, Molecule &mymolecule, bool doprint, std::o
scf_algorithm,scf_alpha,diis_histories,
mygrid->ispin,mygrid->n2ft3d,mymolecule.rho1);


while ((icount < it_out*it_in) && (!converged))
while ((icount < (it_out*it_in)) && (!converged))
{
++icount;
if (stalled)
Expand Down Expand Up @@ -407,7 +406,7 @@ double cgsd_energy(Control2 &control, Molecule &mymolecule, bool doprint, std::o

total_energy = cgsd_bybminimize2(mymolecule,mygeodesic12.mygeodesic1,E,&deltae,
&deltac,bfgscount,ks_it_in,ks_it_out,
scfmix,tole,tolc);
scfmix,tole,tolc,coutput);
++bfgscount;
if (oprint)
coutput << Ifmt(10) << icount
Expand All @@ -416,9 +415,10 @@ double cgsd_energy(Control2 &control, Molecule &mymolecule, bool doprint, std::o
<< Efmt(16,6) << deltac << std::endl;
if ((std::fabs(deltae) > fabs(deltae_old)) ||
(std::fabs(deltae) > 1.0e-2) || (deltae > 0.0))
stalled = true;
stalled = true;
else
stalled = false;
stalled = false;

converged = (std::fabs(deltae) < tole) && (deltac < tolc);
}
} else if (minimizer == 9) {
Expand Down

0 comments on commit 94e7146

Please sign in to comment.