diff --git a/Nwpw/band/CMakeLists.txt b/Nwpw/band/CMakeLists.txt index 31b65c7c..196dfa45 100644 --- a/Nwpw/band/CMakeLists.txt +++ b/Nwpw/band/CMakeLists.txt @@ -2,16 +2,16 @@ # files in the band library file(GLOB_RECURSE src_band_minimizer minimizer/*.hpp minimizer/*.cpp) file(GLOB_RECURSE src_band_cpsd cpsd/*.hpp cpsd/*.cpp) +file(GLOB_RECURSE src_solid lib/solid/*.hpp lib/solid/*.cpp) file(GLOB_RECURSE src_cKinetic lib/cKinetic/*.hpp lib/cKinetic/*.cpp) file(GLOB_RECURSE src_cCoulomb lib/cCoulomb/*.hpp lib/cCoulomb/*.cpp) file(GLOB_RECURSE src_cExchange-Correlation lib/cExchange-Correlation/*.hpp lib/cExchange-Correlation/*.cpp) file(GLOB_RECURSE src_cpsp lib/cpsp/*.hpp lib/cpsp/*.cpp) file(GLOB_RECURSE src_cElectron lib/cElectron/*.hpp lib/cElectron/*.cpp) file(GLOB_RECURSE src_cpsi lib/cpsi/*.hpp lib/cpsi/*.cpp) -file(GLOB_RECURSE src_solid lib/solid/*.hpp lib/solid/*.cpp) # create to the band library -add_library(band ${src_band_minimizer} ${src_band_cpsd} ${src_cKinetic} ${src_cCoulomb} ${src_cExchange-Correlation} ${src_cpsp} ${src_cElectron} ${src_cpsi} ${src_solid}) +add_library(band ${src_band_minimizer} ${src_band_cpsd} ${src_solid} ${src_cKinetic} ${src_cCoulomb} ${src_cExchange-Correlation} ${src_cpsp} ${src_cElectron} ${src_cpsi} ) # interface nwpwlib library to pspw library target_link_libraries(band nwpwlib) diff --git a/Nwpw/band/lib/cElectron/cElectron.cpp b/Nwpw/band/lib/cElectron/cElectron.cpp index 545cce6b..e170886e 100644 --- a/Nwpw/band/lib/cElectron/cElectron.cpp +++ b/Nwpw/band/lib/cElectron/cElectron.cpp @@ -35,8 +35,8 @@ cElectron_Operators::cElectron_Operators(Cneb *mygrid0, cKinetic_Operator *myke0 neall = mygrid->neq[0] + mygrid->neq[1]; /* allocate memory */ - Hpsi = mygrid->g_allocate(1); - psi_r = mygrid->h_allocate(); + Hpsi = mygrid->g_allocate_nbrillq_all(); + psi_r = mygrid->h_allocate_nbrillq_all(); xcp = mygrid->r_nalloc(ispin); xce = mygrid->r_nalloc(ispin); @@ -48,7 +48,7 @@ cElectron_Operators::cElectron_Operators(Cneb *mygrid0, cKinetic_Operator *myke0 vc = mygrid->c_pack_allocate(0); vcall = mygrid->c_pack_allocate(0); - hmltmp = mygrid->m_allocate(-1, 1); + hmltmp = mygrid->w_allocate_nbrillq_all(); omega = mygrid->lattice->omega(); scal1 = 1.0/((double)((mygrid->nx)*(mygrid->ny)*(mygrid->nz))); diff --git a/Nwpw/band/lib/solid/Solid.cpp b/Nwpw/band/lib/solid/Solid.cpp index d0d6b013..a9e3696e 100644 --- a/Nwpw/band/lib/solid/Solid.cpp +++ b/Nwpw/band/lib/solid/Solid.cpp @@ -14,7 +14,7 @@ namespace pwdft { /******************************************** * * - * Molecule::Molecule * + * Solid::Solid * * * ********************************************/ Solid::Solid(char *infilename, bool wvfnc_initialize, Cneb *mygrid0, @@ -28,6 +28,8 @@ Solid::Solid(char *infilename, bool wvfnc_initialize, Cneb *mygrid0, myelectron = myelectron0; mypsp = mypsp0; + nbrillouin = mygrid->nbrillouin; + nbrillq = mygrid->nbrillq; ispin = mygrid->ispin; neall = mygrid->neq[0] + mygrid->neq[1]; ne[0] = mygrid->ne[0]; @@ -38,8 +40,10 @@ Solid::Solid(char *infilename, bool wvfnc_initialize, Cneb *mygrid0, for (int i = 0; i < 60; ++i) E[i] = 0.0; - psi1 = mygrid->g_allocate(1); - psi2 = mygrid->g_allocate(1); + //psi1 = mygrid->g_allocate(1); + //psi2 = mygrid->g_allocate(1); + psi1 = mygrid->g_allocate_nbrillq_all(); + psi2 = mygrid->g_allocate_nbrillq_all(); rho1 = mygrid->r_nalloc(ispin); rho2 = mygrid->r_nalloc(ispin); rho1_all = mygrid->r_nalloc(ispin); @@ -47,9 +51,11 @@ Solid::Solid(char *infilename, bool wvfnc_initialize, Cneb *mygrid0, dng1 = mygrid->c_pack_allocate(0); dng2 = mygrid->c_pack_allocate(0); - lmbda = mygrid->m_allocate(-1, 1); - hml = mygrid->m_allocate(-1, 1); - eig = new double[mygrid->ne[0] + mygrid->ne[1]]; + //lmbda = mygrid->m_allocate(-1, 1); + //hml = mygrid->m_allocate(-1, 1); + hml = mygrid->w_allocate_nbrillq_all(); + lmbda = mygrid->w_allocate_nbrillq_all(); + eig = new double[nbrillq*(ne[0]+ne[1])]; omega = mygrid->lattice->omega(); scal1 = 1.0 / ((double)((mygrid->nx) * (mygrid->ny) * (mygrid->nz))); @@ -57,17 +63,18 @@ Solid::Solid(char *infilename, bool wvfnc_initialize, Cneb *mygrid0, dv = omega * scal1; n2ft3d = (mygrid->n2ft3d); - shift1 = 2 * (mygrid->npack(1)); + //shift1 = 2 * (mygrid->npack(1)); + shift1 = 2 * (mygrid->npack1_max()); shift2 = (mygrid->n2ft3d); - //newpsi = psi_read(mygrid, infilename, wvfnc_initialize, psi1, coutput); + newpsi = cpsi_read(mygrid, infilename, wvfnc_initialize, psi1, coutput); myelectron->gen_vl_potential(); /*---------------------- testing Electron Operators ---------------------- */ - /* + /* double sum1; - std::cout << "Molecule: Testing Electron Operators " << std::endl; + std::cout << "Solid: Testing Electron Operators " << std::endl; myelectron->run(psi1,rho1,dng1,rho1_all); sum1 = myelectron->energy(psi1,rho1,dng1,rho1_all); @@ -77,12 +84,11 @@ Solid::Solid(char *infilename, bool wvfnc_initialize, Cneb *mygrid0, std::cout << "electron ehartr = " << myelectron->ehartree(dng1) << std::endl; std::cout << "electron exc = " << myelectron->exc(rho1_all) << std::endl; std::cout << "electron pxc = " << myelectron->pxc(rho1) - << std::endl; std::cout << "B electron ehartr = " << - myelectron->ehartree(dng1) << std::endl; std::cout << "electron vl_ave = - " << myelectron->vl_ave(dng1) << std::endl; std::cout << "electron vnl ave - = " << myelectron->vnl_ave(psi1) << std::endl; - - */ + << std::endl; + std::cout << "B electron ehartr = " << myelectron->ehartree(dng1) << std::endl; + std::cout << "electron vl_ave = " << myelectron->vl_ave(dng1) << std::endl; + std::cout << "electron vnl ave =" << myelectron->vnl_ave(psi1) << std::endl; + */ /*---------------------- testing Electron Operators ---------------------- */ } diff --git a/Nwpw/band/lib/solid/Solid.hpp b/Nwpw/band/lib/solid/Solid.hpp index c867af73..0cdc1743 100644 --- a/Nwpw/band/lib/solid/Solid.hpp +++ b/Nwpw/band/lib/solid/Solid.hpp @@ -55,7 +55,7 @@ class Solid { double omega, scal2, scal1, dv; - int ispin, ne[2], neall, n2ft3d, shift1, shift2; + int ispin, ne[2], neall, nbrillq, nbrillouin, n2ft3d, shift1, shift2; int nfft[3]; int version = 3; @@ -95,13 +95,13 @@ class Solid { delete[] eig; } - /* write psi molecule */ + /* write psi solid */ void writepsi(char *output_filename, std::ostream &coutput) { //psi_write(mygrid,&version,nfft,mygrid->lattice->unita_ptr(),&ispin,ne, // psi1,output_filename,coutput); } - /* molecule energy */ + /* solid energy */ double energy() { myelectron->run(psi1, rho1, dng1, rho1_all); E[0] = (myelectron->energy(psi1, rho1, dng1, rho1_all) + myewald->energy()); @@ -114,20 +114,22 @@ class Solid { return E[0]; } - /* molecule energy and eigenvalues */ + /* solid energy and eigenvalues */ double energy_eigenvalues() { myelectron->run(psi1, rho1, dng1, rho1_all); E[0] = (myelectron->energy(psi1, rho1, dng1, rho1_all) + myewald->energy()); /* generate eigenvalues */ myelectron->gen_hml(psi1, hml); - mygrid->m_diagonalize(hml, eig); + //mygrid->m_diagonalize(hml, eig); + mygrid->w_diagonalize(hml, eig); return E[0]; } - /* molecule energy and eigenvalues and other energies and en */ - double gen_all_energies() { + /* solid energy and eigenvalues and other energies and en */ + double gen_all_energies() + { myelectron->run(psi1, rho1, dng1, rho1_all); myelectron->gen_energies_en(psi1, rho1, dng1, rho1_all, E, en); @@ -148,7 +150,8 @@ class Solid { /* generate eigenvalues */ myelectron->gen_hml(psi1, hml); - mygrid->m_diagonalize(hml, eig); + //mygrid->m_diagonalize(hml, eig); + mygrid->w_diagonalize(hml, eig); /* generate dipole */ //mypsp->mydipole->gen_dipole(rho1); @@ -171,13 +174,13 @@ class Solid { return ee; } - /* molecule - generate current hamiltonian */ + /* solid - generate current hamiltonian */ void gen_hml() { myelectron->gen_hml(psi1, hml); } - /* molecule - diagonalize the current hamiltonian */ - void diagonalize() { mygrid->m_diagonalize(hml, eig); } + /* solid - diagonalize the current hamiltonian */ + void diagonalize() { mygrid->w_diagonalize(hml, eig); } - /* molecule - call phafacs and gen_vl_potential and semicore */ + /* solid - call phafacs and gen_vl_potential and semicore */ void phafacs_vl_potential_semicore() { mystrfac->phafac(); myewald->phafac(); @@ -298,7 +301,7 @@ class Solid { std::ios init(NULL); init.copyfmt(os); std::string eoln = "\n"; - os << " ============= energy results (Molecule object) =============" << eoln; + os << " ============= energy results (Solid object) =============" << eoln; os << eoln << eoln; os << std::fixed << " number of electrons: spin up= " << std::setw(11) @@ -338,19 +341,35 @@ class Solid { os << " spring bondings : " << Efmt(19,10) << mysolid.E[71] << " (" << Efmt(15,5) << mysolid.E[71]/mysolid.myion->nion << " /ion)" << std::endl; } - os << eoln; - os << eoln; - os << " orbital energy:" << eoln; - int nn = mysolid.ne[0] - mysolid.ne[1]; - double ev = 27.2116; - for (int i=0; iktoindex(nb); + int pk = mysolid.mygrid->ktop(nb); + + int n = mysolid.ne[0] + mysolid.ne[1]; + double tmpeig[n]; + std:memset(tmpeig,0,n*sizeof(double)); + if (pk==mysolid.mygrid->c3db::parall->taskid_k()) + std::memcpy(tmpeig,mysolid.eig+nbq*n,n*sizeof(double)); + mysolid.mygrid->c3db::parall->Vector_SumAll(3,n,tmpeig); + //std::memcpy(tmpeig,mysolid.eig+nbq*n,n*sizeof(double)); + + os << eoln; + os << eoln; + os << mysolid.mygrid->mybrillouin->print_zone_point(nb); + os << eoln; + os << " orbital energies:" << eoln; + int nn = mysolid.ne[0] - mysolid.ne[1]; + double ev = 27.2116; + for (int i=0; imydipole->shortprint_dipole(); diff --git a/Nwpw/band/lib/solid/band_Geodesic.hpp b/Nwpw/band/lib/solid/band_Geodesic.hpp new file mode 100644 index 00000000..6e9cfdda --- /dev/null +++ b/Nwpw/band/lib/solid/band_Geodesic.hpp @@ -0,0 +1,158 @@ +#ifndef _BAND_GEODESIC_HPP_ +#define _BAND_GEODESIC_HPP_ + +#pragma once + +#include "cElectron.hpp" +#include "Solid.hpp" +#include "Cneb.hpp" +#include + +//#include "util.hpp" + +namespace pwdft { + +class band_Geodesic { + + int minimizer; + Solid *mysolid; + cElectron_Operators *myelectron; + + double *U, *Vt, *S; + + // tmp space for ffm multiplies + double *tmp1, *tmp2, *tmp3, *tmpC, *tmpS; + +public: + Cneb *mygrid; + + /* Constructors */ + band_Geodesic(int minimizer0, Solid *mysolid0) { + mysolid = mysolid0; + minimizer = minimizer0; + myelectron = mysolid->myelectron; + mygrid = mysolid->mygrid; + U = mygrid->g_allocate(1); + Vt = mygrid->m_allocate(-1, 1); + S = new double[mygrid->ne[0] + mygrid->ne[1]]; + + // tmp space + tmp1 = mygrid->m_allocate(-1, 1); + tmp2 = mygrid->m_allocate(-1, 1); + tmp3 = mygrid->m_allocate(-1, 1); + tmpC = new double[mygrid->ne[0] + mygrid->ne[1]]; + tmpS = new double[mygrid->ne[0] + mygrid->ne[1]]; + } + + /* destructor */ + ~band_Geodesic() { + delete[] tmpS; + delete[] tmpC; + delete[] tmp3; + delete[] tmp2; + delete[] tmp1; + delete[] S; + delete[] Vt; + delete[] U; + } + + double start(double *A, double *max_sigma, double *min_sigma) { + double *V = mygrid->m_allocate(-1, 1); + //mygrid->ggm_SVD(A, U, S, V); + + int neall = mygrid->ne[0] + mygrid->ne[1]; + double mmsig = 9.99e9; + double msig = 0.0; + for (int i = 0; i < neall; ++i) { + if (std::fabs(S[i]) > msig) + msig = fabs(S[i]); + if (std::fabs(S[i]) < mmsig) + mmsig = fabs(S[i]); + } + *max_sigma = msig; + *min_sigma = mmsig; + + /* calculate Vt */ + mygrid->mm_transpose(-1, V, Vt); + + // double *tmp1 = mygrid->m_allocate(-1,1); + // mygrid->mmm_Multiply(-1,Vt,V,1.0,tmp1,0.0); + // util_matprint("Vt*V",4,tmp1); + + // mygrid->ggm_sym_Multiply(U,U,tmp1); + // util_matprint("Ut*U",4,tmp1); + // delete [] tmp1; + + delete[] V; + + /* calculate and return 2* */ + return (2.0 * myelectron->eorbit(A)); + } + + + void get(double t, double *Yold, double *Ynew) { + //mygrid->mm_SCtimesVtrans(-1, t, S, Vt, tmp1, tmp3, tmpC, tmpS); + + /* Ynew = Yold*V*cos(Sigma*t)*Vt + U*sin(Sigma*t)*Vt */ + mygrid->mmm_Multiply2(-1, Vt, tmp1, 1.0, tmp2, 0.0); + //mygrid->fmf_Multiply(-1, Yold, tmp2, 1.0, Ynew, 0.0); + //mygrid->fmf_Multiply(-1, U, tmp3, 1.0, Ynew, 1.0); + + /* ortho check - need to figure out what causes this to happen */ + double sum2 = mygrid->gg_traceall(Ynew, Ynew); + double sum1 = mygrid->ne[0] + mygrid->ne[1]; + if ((mygrid->ispin) == 1) + sum1 *= 2; + if (std::fabs(sum2 - sum1) > 1.0e-10) { + // if (myparall->is_master()) std::cout << " Warning - Gram-Schmidt being + // performed on psi2" << std::endl; std::cout << " Warning - Gram-Schmidt + // being performed on psi2, t=" + // << t << " sum1=" << sum1 << " sum2=" << sum2 << " error=" << + // fabs(sum2-sum1) << std::endl; + mygrid->g_ortho(Ynew); + } + } + + void transport(double t, double *Yold, double *Ynew) { + //mygrid->mm_SCtimesVtrans2(-1, t, S, Vt, tmp1, tmp3, tmpC, tmpS); + + /* tHnew = (-Yold*V*sin(Sigma*t) + U*cos(Sigma*t))*Sigma*Vt */ + mygrid->mmm_Multiply2(-1, Vt, tmp1, 1.0, tmp2, 0.0); + //mygrid->fmf_Multiply(-1, Yold, tmp2, -1.0, Ynew, 0.0); + //mygrid->fmf_Multiply(-1, U, tmp3, 1.0, Ynew, 1.0); + } + + void psi_1transport(double t, double *H0) { + this->transport(t, mysolid->psi1, H0); + } + + void Gtransport(double t, double *Yold, double *tG) { + // mygrid->ffm_sym_Multiply(-1,U,tG,tmp2); + //mygrid->ffm_Multiply(-1, U, tG, tmp2); + //mygrid->mm_SCtimesVtrans3(-1, t, S, tmp2, tmp1, tmp3, tmpC, tmpS); + mygrid->mmm_Multiply2(-1, Vt, tmp1, 1.0, tmp2, 0.0); + + //mygrid->fmf_Multiply(-1, Yold, tmp2, -1.0, tG, 1.0); + //mygrid->fmf_Multiply(-1, U, tmp3, -1.0, tG, 1.0); + } + + void psi_1Gtransport(double t, double *H0) { + this->Gtransport(t, mysolid->psi1, H0); + } + + double energy(double t) { + this->get(t, mysolid->psi1, mysolid->psi2); + return (mysolid->psi2_energy()); + } + + double denergy(double t) { + this->transport(t, mysolid->psi1, mysolid->psi2); + return (2.0 * mysolid->psi2_eorbit()); + } + + void psi_final(double t) { this->get(t, mysolid->psi1, mysolid->psi2); } +}; + +} // namespace pwdft + +#endif diff --git a/Nwpw/band/lib/solid/band_Geodesic12.cpp b/Nwpw/band/lib/solid/band_Geodesic12.cpp new file mode 100644 index 00000000..522aded4 --- /dev/null +++ b/Nwpw/band/lib/solid/band_Geodesic12.cpp @@ -0,0 +1,27 @@ +/* band_Geodesic12.cpp - + Author - Eric Bylaska +*/ + +#include "band_Geodesic12.hpp" + +namespace pwdft { + +/******************************************* + * * + * band_Geodesic12::band_Geodesic12 * + * * + *******************************************/ +band_Geodesic12::band_Geodesic12(int minimizer0, Solid *mysolid0, Control2 &control) { + int minimizer = control.minimizer(); + if ((minimizer == 1) || (minimizer == 2)) { + has_geodesic1 = true; + mygeodesic1 = new (std::nothrow) band_Geodesic(minimizer0, mysolid0); + } + + if ((minimizer == 4) || (minimizer == 7)) { + has_geodesic2 = true; + mygeodesic2 = new (std::nothrow) band_Geodesic2(minimizer0, mysolid0); + } +} + +} // namespace pwdft diff --git a/Nwpw/band/lib/solid/band_Geodesic12.hpp b/Nwpw/band/lib/solid/band_Geodesic12.hpp new file mode 100644 index 00000000..33e8ba56 --- /dev/null +++ b/Nwpw/band/lib/solid/band_Geodesic12.hpp @@ -0,0 +1,31 @@ +#pragma once + +#include "Control2.hpp" +#include "band_Geodesic.hpp" +#include "band_Geodesic2.hpp" +#include "Cneb.hpp" + +namespace pwdft { + +class band_Geodesic12 { + +public: + bool has_geodesic1 = false; + band_Geodesic *mygeodesic1; + + bool has_geodesic2 = false; + band_Geodesic2 *mygeodesic2; + + /* Constructors */ + band_Geodesic12(int, Solid *, Control2 &); + + /* destructor */ + ~band_Geodesic12() { + if (has_geodesic1) + delete mygeodesic1; + if (has_geodesic2) + delete mygeodesic2; + } +}; + +} // namespace pwdft diff --git a/Nwpw/band/lib/solid/band_Geodesic2.hpp b/Nwpw/band/lib/solid/band_Geodesic2.hpp new file mode 100644 index 00000000..eceb3f15 --- /dev/null +++ b/Nwpw/band/lib/solid/band_Geodesic2.hpp @@ -0,0 +1,215 @@ +#ifndef _BAND_GEODESIC2_HPP_ +#define _BAND_GEODESIC2_HPP_ + +#pragma once + +#include "cElectron.hpp" +#include "Solid.hpp" +#include "Cneb.hpp" +#include + +//#include "util.hpp" + +namespace pwdft { + +class band_Geodesic2 { + + int minimizer; + Solid *mysolid; + cElectron_Operators *myelectron; + + // tmp space of m matrices + double *R, *A, *MM, *NN, *TT; + + // tmp space of m4 matrices + double *T4, *V4, *W4, *A4, *B4, *R4; + + // tmp space for ffm multiplies + double *Q, *Hold, *S; + int nemax, npack1; + +public: + Cneb *mygrid; + + /* Constructors */ + band_Geodesic2(int minimizer0, Solid *mysolid0) { + mysolid = mysolid0; + minimizer = minimizer0; + myelectron = mysolid->myelectron; + mygrid = mysolid->mygrid; + + nemax = mygrid->neq[0] + mygrid->neq[1]; + npack1 = mygrid->CGrid::npack(1); + + Hold = mygrid->g_allocate(1); + Q = mygrid->g_allocate(1); + + R = mygrid->m_allocate(-1, 1); + A = mygrid->m_allocate(-1, 1); + MM = mygrid->m_allocate(-1, 1); + NN = mygrid->m_allocate(-1, 1); + TT = mygrid->m_allocate(-1, 1); + + //T4 = mygrid->m4_allocate(-1, 1); + //V4 = mygrid->m4_allocate(-1, 1); + //W4 = mygrid->m4_allocate(-1, 1); + //A4 = mygrid->m4_allocate(-1, 1); + //B4 = mygrid->m4_allocate(-1, 1); + //R4 = mygrid->m4_allocate(-1, 1); + + S = new double[2 * (mygrid->ne[0] + mygrid->ne[1])]; + // tmp space + } + + /* destructor */ + ~band_Geodesic2() { + delete[] Hold; + delete[] Q; + + delete[] R; + delete[] A; + delete[] MM; + delete[] NN; + + delete[] T4; + delete[] V4; + delete[] W4; + delete[] A4; + delete[] B4; + delete[] R4; + + delete[] S; + } + + /***************************************** + * * + * start * + * * + ***************************************** + + This function determines the pxp matrices R and YA, and + the orthogonal nxp matrix Q. Q and R are determined from + the QR decomposition of the projected direction (I-YY^t)H, and + YH is defined as the Lagrange Multiplier pxp matrix Y^tH. + */ + + double start(double *Y, double *H, double *max_sigma, double *min_sigma) { + + int nn = 2 * npack1 * nemax; + int one = 1; + double mrone = 1.0; + + //**** Hold <-- H **** + std::memcpy(Hold, H, 2 * nemax * npack1 * sizeof(double)); + + //**** calculate A= **** + //mygrid->ffm_Multiply(-1, Y, H, A); + + //**** calculate Q=(I-YYt)H - should not be necessary but just in case **** + //mygrid->fmf_Multiply(-1, Y, A, 1.0, Q, 0.0); + DAXPY_PWDFT(nn, mrone, H, one, Q, one); + DSCAL_PWDFT(nn, mrone, Q, one); + + //**** calculate QR using Modified Gram-Schmidt **** + mygrid->fm_QR(-1, Q, R); + + //**** generate T from A and R **** + //mygrid->mmm4_AR_to_T4(-1, A, R, T4); + + //**** Factor T--> V,W,and S **** + // mygrid->m4_FactorSkew(-1, T4, V4, W4, S); + + /* calculate and return 2* */ + return (2.0 * myelectron->eorbit(H)); + } + + /***************************************** + * * + * get * + * * + ***************************************** + This routine calculates + + Ynew = Yold*M(t) + Q*N(t) + + where + - - - - + | M(t) | = Exp(t*T)* | I | + | N(t) | | 0 | + - - - - + */ + void get(double t, double *Yold, double *Ynew) { + this->get_MandN(t, MM, NN); + //mygrid->fmf_Multiply(-1, Yold, MM, 1.0, Ynew, 0.0); + //mygrid->fmf_Multiply(-1, Q, NN, 1.0, Ynew, 1.0); + } + + /***************************************** + * * + * transport * + * * + ***************************************** + This routine calculates + + Hnew = Hold*M(t) + Yold*R^t*N(t) + + where + - - - - + | M(t) | = Exp(t*T)* | I | + | N(t) | | 0 | + - - - - + */ + void transport(double t, double *Yold, double *Hnew) { + + this->get_MandN(t, MM, NN); + + //*** TT(t) = -R^t*NN(t) *** + //mygrid->mmm_Multiply2(-1, R, NN, -1.0, TT, 0.0); + + //*** Hnew <-- Hold*M(t) + Yold*TT(t) *** + //mygrid->fmf_Multiply(-1, Hold, MM, 1.0, Hnew, 0.0); + //mygrid->fmf_Multiply(-1, Yold, TT, 1.0, Hnew, 1.0); + } + + /***************************************** + * * + * get_MandN * + * * + ***************************************** + This routine returns + - - - - + | M(t) | = Exp(t*T)* | I | + | N(t) | | 0 | + - - - - + where + + T = U*Sigma*U^H, with U=(V+iW) + + is a skew matrix that is decomposed into V,W,and Sigma + + */ + void get_MandN(const double t, double *M, double *N) { + //mygrid->m4_RotationSkew(-1, t, V4, W4, S, A4, B4, R4); + //mygrid->m4_R4_to_MN(-1, R4, M, N); + } + + void psi_1transport(double t, double *H0) { + this->transport(t, mysolid->psi1, H0); + } + + double energy(double t) { + this->get(t, mysolid->psi1, mysolid->psi2); + return (mysolid->psi2_energy()); + } + + double denergy(double t) { + this->transport(t, mysolid->psi1, mysolid->psi2); + return (2.0 * mysolid->psi2_eorbit()); + } + + void psi_final(double t) { this->get(t, mysolid->psi1, mysolid->psi2); } +}; + +} // namespace pwdft + +#endif diff --git a/Nwpw/band/lib/solid/band_lmbfgs.hpp b/Nwpw/band/lib/solid/band_lmbfgs.hpp new file mode 100644 index 00000000..cce10431 --- /dev/null +++ b/Nwpw/band/lib/solid/band_lmbfgs.hpp @@ -0,0 +1,134 @@ +#ifndef _BAND_LMBFGS_HPP_ +#define _BAND_LMBFGS_HPP_ + +#pragma once + +#include "band_Geodesic.hpp" +#include + +namespace pwdft { + +class band_lmbfgs { + + band_Geodesic *mygeodesic; + + int npack1, neall, nsize, max_m, m; + int indx[20], size_list; + double rho[20]; + double *lm_list; + +public: + /* Constructors */ + band_lmbfgs(band_Geodesic *mygeodesic0, const int max_m0) { + + mygeodesic = mygeodesic0; + max_m = max_m0; + npack1 = mygeodesic->mygrid->npack(1); + neall = mygeodesic->mygrid->neq[0] + mygeodesic->mygrid->neq[1]; + nsize = 2 * neall * npack1; + size_list = 2 * max_m; + for (auto k = 0; k < size_list; ++k) + indx[k] = k; + + m = 0; + + lm_list = mygeodesic->mygrid->g_nallocate(1, size_list); + + // std::memcpy(&lm_list[2*m*nsize],g,nsize*sizeof(double)); + // mygeodesic->mygrid->g_Scale(-1.0,&lm_list[2*m*nsize]); + // std::memcpy(&lm_list[(2*m+1)*nsize],&lm_list[2*m*nsize],nsize*sizeof(double)); + } + + /* destructor */ + ~band_lmbfgs() { mygeodesic->mygrid->g_deallocate(lm_list); } + + void start(const double *g) { + std::memcpy(&lm_list[2 * m * nsize], g, nsize * sizeof(double)); + mygeodesic->mygrid->g_Scale(-1.0, &lm_list[2 * m * nsize]); + std::memcpy(&lm_list[(2 * m + 1) * nsize], &lm_list[2 * m * nsize], + nsize * sizeof(double)); + } + + void fetch(const double tmin, double *g, double *s) { + double *yy, *ss, sum, sum0, alpha[20], beta; + + mygeodesic->mygrid->g_Scale(-1.0, g); + + std::memcpy(s, g, nsize * sizeof(double)); + + yy = &lm_list[indx[2 * m] * nsize]; + ss = &lm_list[indx[2 * m + 1] * nsize]; + + mygeodesic->psi_1Gtransport(tmin, yy); + mygeodesic->psi_1transport(tmin, ss); + + mygeodesic->mygrid->gg_daxpy(-1.0, g, yy); + mygeodesic->mygrid->g_Scale(-1.0, yy); + + sum = mygeodesic->mygrid->gg_traceall(yy, ss); + + //*** exit if dividing by small number *** + if (std::abs(sum) > 1.0e-15) { + rho[m] = 1.0 / sum; + + sum = mygeodesic->mygrid->gg_traceall(ss, s); + alpha[m] = rho[m] * sum; + + mygeodesic->mygrid->gg_daxpy(-alpha[m], yy, s); + for (auto k = m - 1; k > 0; --k) { + yy = &lm_list[indx[2 * (k + 1)] * nsize]; + ss = &lm_list[indx[2 * (k + 1) + 1] * nsize]; + + mygeodesic->psi_1Gtransport(tmin, yy); + mygeodesic->psi_1Gtransport(tmin, ss); + + sum = mygeodesic->mygrid->gg_traceall(ss, s); + alpha[k] = rho[k] * sum; + + mygeodesic->mygrid->gg_daxpy(-alpha[k], yy, s); + } + + //**** preconditioner **** + // call Grsm_gg_dScale(npack1,neall,h0,s,s) + + for (auto k = 0; k < (m - 1); ++k) { + sum = mygeodesic->mygrid->gg_traceall(yy, s); + beta = rho[k] * sum; + sum0 = alpha[k] - beta; + + mygeodesic->mygrid->gg_daxpy(sum0, ss, s); + + yy = &lm_list[indx[2 * (k + 1)] * nsize]; + ss = &lm_list[indx[2 * (k + 1) + 1] * nsize]; + } + + sum = mygeodesic->mygrid->gg_traceall(yy, s); + beta = rho[m] * sum; + sum0 = alpha[m] - beta; + + mygeodesic->mygrid->gg_daxpy(sum0, ss, s); + + if (m < (max_m - 1)) + ++m; + else { + int itmp0 = indx[0]; + int itmp1 = indx[1]; + for (auto k = 0; k < (size_list - 2); ++k) + indx[k] = indx[k + 2]; + indx[size_list - 2] = itmp0; + indx[size_list - 1] = itmp1; + + for (auto k = 0; k < (m - 1); ++k) + rho[k] = rho[k + 1]; + } + mygeodesic->mygrid->g_Scale(-1.0, s); + } + std::memcpy(&lm_list[indx[2 * m] * nsize], g, nsize * sizeof(double)); + std::memcpy(&lm_list[indx[2 * m + 1] * nsize], s, nsize * sizeof(double)); + mygeodesic->mygrid->g_Scale(-1.0, g); + } +}; + +} // namespace pwdft + +#endif diff --git a/Nwpw/band/lib/solid/band_lmbfgs2.hpp b/Nwpw/band/lib/solid/band_lmbfgs2.hpp new file mode 100644 index 00000000..3b5d3ffa --- /dev/null +++ b/Nwpw/band/lib/solid/band_lmbfgs2.hpp @@ -0,0 +1,134 @@ +#ifndef _BAND_LMBFGS2_HPP_ +#define _BAND_LMBFGS2_HPP_ + +#pragma once + +#include "band_Geodesic2.hpp" +#include + +namespace pwdft { + +class band_lmbfgs2 { + + band_Geodesic2 *mygeodesic; + + int npack1, neall, nsize, max_m, m; + int indx[20], size_list; + double rho[20]; + double *lm_list; + +public: + /* Constructors */ + band_lmbfgs2(band_Geodesic2 *mygeodesic0, const int max_m0) { + + mygeodesic = mygeodesic0; + max_m = max_m0; + npack1 = mygeodesic->mygrid->npack(1); + neall = mygeodesic->mygrid->neq[0] + mygeodesic->mygrid->neq[1]; + nsize = 2 * neall * npack1; + size_list = 2 * max_m; + for (auto k = 0; k < size_list; ++k) + indx[k] = k; + + m = 0; + + lm_list = mygeodesic->mygrid->g_nallocate(1, size_list); + + // std::memcpy(&lm_list[2*m*nsize],g,nsize*sizeof(double)); + // mygeodesic->mygrid->g_Scale(-1.0,&lm_list[2*m*nsize]); + // std::memcpy(&lm_list[(2*m+1)*nsize],&lm_list[2*m*nsize],nsize*sizeof(double)); + } + + /* destructor */ + ~band_lmbfgs2() { mygeodesic->mygrid->g_deallocate(lm_list); } + + void start(const double *g) { + std::memcpy(&lm_list[2 * m * nsize], g, nsize * sizeof(double)); + mygeodesic->mygrid->g_Scale(-1.0, &lm_list[2 * m * nsize]); + std::memcpy(&lm_list[(2 * m + 1) * nsize], &lm_list[2 * m * nsize], + nsize * sizeof(double)); + } + + void fetch(const double tmin, double *g, double *s) { + double *yy, *ss, sum, sum0, alpha[20], beta; + + mygeodesic->mygrid->g_Scale(-1.0, g); + + std::memcpy(s, g, nsize * sizeof(double)); + + yy = &lm_list[indx[2 * m] * nsize]; + ss = &lm_list[indx[2 * m + 1] * nsize]; + + // mygeodesic->psi_1Gtransport(tmin,yy); + mygeodesic->psi_1transport(tmin, ss); + + mygeodesic->mygrid->gg_daxpy(-1.0, g, yy); + mygeodesic->mygrid->g_Scale(-1.0, yy); + + sum = mygeodesic->mygrid->gg_traceall(yy, ss); + + //*** exit if dividing by small number *** + if (std::abs(sum) > 1.0e-15) { + rho[m] = 1.0 / sum; + + sum = mygeodesic->mygrid->gg_traceall(ss, s); + alpha[m] = rho[m] * sum; + + mygeodesic->mygrid->gg_daxpy(-alpha[m], yy, s); + for (auto k = m - 1; k > 0; --k) { + yy = &lm_list[indx[2 * (k + 1)] * nsize]; + ss = &lm_list[indx[2 * (k + 1) + 1] * nsize]; + + // mygeodesic->psi_1Gtransport(tmin,yy); + // mygeodesic->psi_1Gtransport(tmin,ss); + + sum = mygeodesic->mygrid->gg_traceall(ss, s); + alpha[k] = rho[k] * sum; + + mygeodesic->mygrid->gg_daxpy(-alpha[k], yy, s); + } + + //**** preconditioner **** + // call Grsm_gg_dScale(npack1,neall,h0,s,s) + + for (auto k = 0; k < (m - 1); ++k) { + sum = mygeodesic->mygrid->gg_traceall(yy, s); + beta = rho[k] * sum; + sum0 = alpha[k] - beta; + + mygeodesic->mygrid->gg_daxpy(sum0, ss, s); + + yy = &lm_list[indx[2 * (k + 1)] * nsize]; + ss = &lm_list[indx[2 * (k + 1) + 1] * nsize]; + } + + sum = mygeodesic->mygrid->gg_traceall(yy, s); + beta = rho[m] * sum; + sum0 = alpha[m] - beta; + + mygeodesic->mygrid->gg_daxpy(sum0, ss, s); + + if (m < (max_m - 1)) + ++m; + else { + int itmp0 = indx[0]; + int itmp1 = indx[1]; + for (auto k = 0; k < (size_list - 2); ++k) + indx[k] = indx[k + 2]; + indx[size_list - 2] = itmp0; + indx[size_list - 1] = itmp1; + + for (auto k = 0; k < (m - 1); ++k) + rho[k] = rho[k + 1]; + } + mygeodesic->mygrid->g_Scale(-1.0, s); + } + std::memcpy(&lm_list[indx[2 * m] * nsize], g, nsize * sizeof(double)); + std::memcpy(&lm_list[indx[2 * m + 1] * nsize], s, nsize * sizeof(double)); + mygeodesic->mygrid->g_Scale(-1.0, g); + } +}; + +} // namespace pwdft + +#endif diff --git a/Nwpw/band/minimizer/band_cgsd.hpp b/Nwpw/band/minimizer/band_cgsd.hpp new file mode 100644 index 00000000..3566d412 --- /dev/null +++ b/Nwpw/band/minimizer/band_cgsd.hpp @@ -0,0 +1,25 @@ +#ifndef _BAND_CGSD_HPP_ +#define _BAND_CGSD_HPP_ + +#pragma once + +namespace pwdft { + +#include "band_Geodesic.hpp" +#include "band_Geodesic2.hpp" +#include "Solid.hpp" +#include "band_lmbfgs.hpp" + +extern double band_cgsd_cgminimize(Solid &, band_Geodesic *, double *, double *, + double *, int, int, double, double); +extern double band_cgsd_bfgsminimize(Solid &, band_Geodesic *, band_lmbfgs &, double *, + double *, double *, int, int, double, double); + +extern double band_cgsd_cgminimize2(Solid &, band_Geodesic2 *, double *, double *, + double *, int, int, double, double); +extern double band_cgsd_bfgsminimize2(Solid &, band_Geodesic2 *, band_lmbfgs2 &, + double *, double *, double *, int, int, double, + double); + +} // namespace pwdft +#endif diff --git a/Nwpw/band/minimizer/band_cgsd_cgminimize.cpp b/Nwpw/band/minimizer/band_cgsd_cgminimize.cpp new file mode 100644 index 00000000..04775e7f --- /dev/null +++ b/Nwpw/band/minimizer/band_cgsd_cgminimize.cpp @@ -0,0 +1,128 @@ + +#include +#include +#include +#include +#include + +#include "Control2.hpp" +#include "band_Geodesic.hpp" +#include "Ion.hpp" +#include "Solid.hpp" +#include "Parallel.hpp" +#include "Cneb.hpp" +#include "util_date.hpp" +#include "util_linesearch.hpp" + +namespace pwdft { + +/* create dummy function call to Geodesic class functions */ +static band_Geodesic *mygeodesic_ptr; +static double dummy_energy(double t) { return mygeodesic_ptr->energy(t); } +static double dummy_denergy(double t) { return mygeodesic_ptr->denergy(t); } + +/****************************************** + * * + * band_cgsd_cgminimize * + * * + ******************************************/ +double band_cgsd_cgminimize(Solid &mysolid, band_Geodesic *mygeodesic, double *E, + double *deltae, double *deltac, int current_iteration, + int it_in, double tole, double tolc) { + bool done = false; + double tmin = 0.0; + double deltat_min = 1.0e-3; + double deltat; + double sum0, sum1, scale, total_energy; + double dE, max_sigma, min_sigma; + double Eold, dEold, Enew; + double tmin0, deltae0; + + Cneb *mygrid = mysolid.mygrid; + mygeodesic_ptr = mygeodesic; + + /* get the initial gradient and direction */ + double *G1 = mygrid->g_allocate(1); + double *H0 = mygrid->g_allocate(1); + + //|-\____|\/-----\/\/-> Start Parallel Section <-\/\/-----\/|____/-| + + total_energy = mysolid.psi_1get_Tgradient(G1); + sum1 = mygrid->gg_traceall(G1, G1); + Enew = total_energy; + + mygrid->gg_copy(G1, H0); + + /****************************************** + **** **** + **** Start of conjugate gradient loop **** + **** **** + ******************************************/ + int it = 0; + tmin = deltat_min; + while ((!done) && ((it++) < it_in)) { + /* initialize the geoedesic line data structure */ + dEold = mygeodesic->start(H0, &max_sigma, &min_sigma); + + /* line search */ + if (tmin > deltat_min) + deltat = tmin; + else + deltat = deltat_min; + + tmin0 = tmin; + deltae0 = *deltae; + + Eold = Enew; + Enew = util_linesearch(0.0, Eold, dEold, deltat, &dummy_energy, + &dummy_denergy, 0.50, &tmin0, &deltae0, 2); + tmin = tmin0; + *deltae = deltae0; + *deltac = mysolid.rho_error(); + mygeodesic->psi_final(tmin); + + /* exit loop early */ + done = ((it >= it_in) || ((std::fabs(*deltae) < tole) && (*deltac < tolc))); + + /* transport the previous search directions */ + mygeodesic->psi_1transport(tmin, H0); + + /* make psi1 <--- psi2(tmin) */ + mysolid.swap_psi1_psi2(); + + if (!done) { + /* get the new gradient - also updates densities */ + total_energy = mysolid.psi_1get_Tgradient(G1); + sum0 = sum1; + sum1 = mygrid->gg_traceall(G1, G1); + + /* the new direction using Fletcher-Reeves */ + if ((std::fabs(*deltae) <= (1.0e-2)) && (tmin > deltat_min)) { + if (sum0 > 1.0e-9) + scale = sum1 / sum0; + else + scale = 0.0; + + mygrid->g_Scale(scale, H0); + mygrid->gg_Sum2(G1, H0); + } + + /* the new direction using steepest-descent */ + else + mygrid->gg_copy(G1, H0); + + // mygrid->gg_copy(G1,H0); + } + } + // Making an extra call to electron.run and energy + total_energy = mysolid.gen_all_energies(); + + //|-\____|\/-----\/\/-> End Parallel Section <-\/\/-----\/|____/-| + + mygrid->g_deallocate(H0); + mygrid->g_deallocate(G1); + + return total_energy; +} + +} // namespace pwdft diff --git a/Nwpw/band/minimizer/band_cgsd_energy.cpp b/Nwpw/band/minimizer/band_cgsd_energy.cpp new file mode 100644 index 00000000..2c80f767 --- /dev/null +++ b/Nwpw/band/minimizer/band_cgsd_energy.cpp @@ -0,0 +1,283 @@ +#include +#include +#include +#include +#include + +#include "Control2.hpp" +#include "band_Geodesic12.hpp" +#include "Ion.hpp" +#include "Solid.hpp" +#include "Parallel.hpp" +#include "Cneb.hpp" +#include "iofmt.hpp" +#include "band_lmbfgs.hpp" +#include "band_lmbfgs2.hpp" +#include "util_date.hpp" + +#include "band_cgsd.hpp" + +#include "iofmt.hpp" + +namespace pwdft { + +/****************************************** + * * + * band_cgsd_noit_energy * + * * + ******************************************/ +double band_cgsd_noit_energy(Solid &mysolid, bool doprint, std::ostream &coutput) +{ + Parallel *parall = mysolid.mygrid->c3db::parall; + + /* generate phase factors and local psp and semicore density */ + mysolid.phafacs_vl_potential_semicore(); + + double total_energy = mysolid.gen_all_energies(); + + /* report summary of results */ + if (parall->base_stdio_print && doprint) + { + coutput << " ================== optimization turned off ===================" << std::endl << std::endl; + coutput << std::endl; + coutput << mysolid; + } + + return total_energy; +} + +/****************************************** + * * + * band_cgsd_energy * + * * + ******************************************/ +double band_cgsd_energy(Control2 &control, Solid &mysolid, bool doprint, std::ostream &coutput) +{ + Parallel *parall = mysolid.mygrid->c3db::parall; + Cneb *mygrid = mysolid.mygrid; + Ion *myion = mysolid.myion; + + bool stalled = false; + int ne[2],ispin,nion; + double E[70],total_energy,deltae,deltae_old,deltac; + + int it_in = control.loop(0); + int it_out = control.loop(1); + double tole = control.tolerances(0); + double tolc = control.tolerances(1); + + double dt = control.time_step(); + double dte = dt/sqrt(control.fake_mass()); + + int minimizer = control.minimizer(); + int lmbfgs_size = control.lmbfgs_size(); + + bool hprint = (parall->is_master() && control.print_level("high") && doprint); + bool oprint = (parall->is_master() && control.print_level("medium") && doprint); + bool lprint = (parall->is_master() && control.print_level("low") && doprint); + + for (auto ii=0; ii<70; ++ii) + E[ii] = 0.0; + + if (oprint) + { + if (minimizer == 1) coutput << " =========== Grassmann conjugate gradient iteration ===========" << std::endl; + if (minimizer == 2) coutput << " ================= Grassmann lmbfgs iteration =================" << std::endl; + if (minimizer == 4) coutput << " ============ Stiefel conjugate gradient iteration ============" << std::endl; + if (minimizer == 5) coutput << " ============ Kohn-Sham scf iteration (potential) =============" << std::endl; + if (minimizer == 7) coutput << " ================== Stiefel lmbfgs iteration ==================" << std::endl; + if (minimizer == 8) coutput << " ============= Kohn-Sham scf iteration (density) ==============" << std::endl; + + coutput << " >>> iteration started at " << util_date() << " <<<" << std::endl;; + coutput << " iter. Energy DeltaE DeltaRho" << std::endl; + coutput << " --------------------------------------------------------------" << std::endl; + // printf("%10d%25.12le%16.6le%16.6le\n",1000,99.99, 1.33434340e-4, 2.33434211e-6); + } + + // if (minimizer > 1) band_Grsm_list_start() + if ((minimizer == 5) || (minimizer == 8)) + it_out = 1; + + band_Geodesic12 mygeodesic12(minimizer, &mysolid, control); + + /* generate phase factors and local psp and semicore density */ + mysolid.phafacs_vl_potential_semicore(); + + // std::cout << "band_cgsd_energy: minimizer = " << minimizer << std::endl; + deltae = -1.0e-03; + int bfgscount = 0; + int icount = 0; + bool converged = false; + + if (minimizer == 1) { + if (mysolid.newpsi) + { + int it_in0 = 15; + for (int it=0; it fabs(deltae_old)) || + (std::fabs(deltae) > 1.0e-2) || (deltae > 0.0)) + stalled = true; + else + stalled = false; + converged = (std::fabs(deltae) < tole) && (deltac < tolc); + } + + } else if (minimizer == 2) { + if (mysolid.newpsi) { + int it_in0 = 15; + for (int it=0; it fabs(deltae_old)) || + (std::fabs(deltae) > 1.0e-2) || (deltae > 0.0)) + stalled = true; + else + stalled = false; + converged = (std::fabs(deltae) < tole) && (deltac < tolc); + } + } else if (minimizer == 4) { + if (mysolid.newpsi) { + int it_in0 = 15; + for (int it = 0; it < it_in0; ++it) + mysolid.sd_update_sic(dte); + if (oprint) + coutput << " - " << it_in0 + << " steepest descent iterations performed" << std::endl; + } + while ((icount < it_out) && (!converged)) { + ++icount; + if (stalled) { + for (int it = 0; it < it_in; ++it) + mysolid.sd_update_sic(dte); + if (oprint) + std::cout << " - " << it_in + << " steepest descent iterations performed" << std::endl; + bfgscount = 0; + } + deltae_old = deltae; + // total_energy = + // band_cgsd_cgminimize2(mysolid, mygeodesic12.mygeodesic2, E, &deltae, + // &deltac, bfgscount, it_in, tole, tolc); + ++bfgscount; + if (oprint) + coutput << Ifmt(10) << icount*it_in + << Efmt(25,12) << total_energy + << Efmt(16,6) << deltae + << Efmt(16,6) << deltac << std::endl; + if ((std::fabs(deltae) > fabs(deltae_old)) || + (std::fabs(deltae) > 1.0e-2) || (deltae > 0.0)) + stalled = true; + else + stalled = false; + converged = (std::fabs(deltae) < tole) && (deltac < tolc); + } + } else if (minimizer == 7) { + if (mysolid.newpsi) { + int it_in0 = 15; + for (int it=0; it fabs(deltae_old)) || + (std::fabs(deltae) > 1.0e-2) || (deltae > 0.0)) + stalled = true; + else + stalled = false; + converged = (std::fabs(deltae) < tole) && (deltac < tolc); + } + } + + if (oprint) { + if (converged) coutput << " *** tolerance ok. iteration terminated" << std::endl; + if (!converged) coutput << " *** arrived at the Maximum iteration. terminated" << std::endl; + coutput << " >>> iteration ended at " << util_date() << " <<<" << std::endl; + } + + /* report summary of results */ + // total_energy = mysolid.gen_all_energies(); + if (oprint) { + coutput << std::endl; + coutput << mysolid; + } + + return total_energy; +} + +/****************************************** + * * + * band_cgsd_energy_gradient * + * * + ******************************************/ +void band_cgsd_energy_gradient(Solid &mysolid, double *grad_ion) +{ + + mysolid.psi_1local_force(grad_ion); + mysolid.psi_1nonlocal_force(grad_ion); + mysolid.semicore_force(grad_ion); + + mysolid.ion_fion(grad_ion); +} + +} // namespace pwdft diff --git a/Nwpw/band/minimizer/band_cgsd_energy.hpp b/Nwpw/band/minimizer/band_cgsd_energy.hpp new file mode 100644 index 00000000..10bc70c4 --- /dev/null +++ b/Nwpw/band/minimizer/band_cgsd_energy.hpp @@ -0,0 +1,15 @@ +#ifndef _BAND_CGSD_ENERGY_HPP_ +#define _BAND_CGSD_ENERGY_HPP_ + +#pragma once + +namespace pwdft { + +#include "Solid.hpp" + +extern double band_cgsd_noit_energy(Solid &, bool, std::ostream &); +extern double band_cgsd_energy(Control2 &, Solid &, bool, std::ostream &); +extern void band_cgsd_energy_gradient(Solid &, double *); + +} // namespace pwdft +#endif diff --git a/Nwpw/band/minimizer/band_minimizer.cpp b/Nwpw/band/minimizer/band_minimizer.cpp index 0d8692ce..48c28123 100644 --- a/Nwpw/band/minimizer/band_minimizer.cpp +++ b/Nwpw/band/minimizer/band_minimizer.cpp @@ -22,6 +22,8 @@ #include "CPseudopotential.hpp" #include "cpsi.hpp" #include "CStrfac.hpp" +#include "Solid.hpp" +#include "cElectron.hpp" #include "util_date.hpp" //#include "rtdb.hpp" @@ -32,7 +34,7 @@ #include "psp_file_check.hpp" #include "psp_library.hpp" -//#include "cgsd_energy.hpp" +#include "band_cgsd_energy.hpp" #include "json.hpp" using json = nlohmann::json; @@ -82,7 +84,7 @@ int band_minimizer(MPI_Comm comm_world0, std::string &rtdbstring, std::ostream & coutput << " * *\n"; coutput << " * PWDFT BAND Calculation *\n"; coutput << " * *\n"; - coutput << " * [ (Grassmann/Stiefel manifold implementation) ] *\n"; + coutput << " * [ (bundled Grassmann/Stiefel manifold) ] *\n"; coutput << " * [ C++ implementation ] *\n"; coutput << " * *\n"; coutput << " * version #7.00 07/19/24 *\n"; @@ -146,6 +148,8 @@ int band_minimizer(MPI_Comm comm_world0, std::string &rtdbstring, std::ostream & // initialize gdevice memory mygrid.c3db::mygdevice.psi_alloc(mygrid.npack1_max(),mygrid.neq[0]+mygrid.neq[1],control.tile_factor()); + + // setup structure factor CStrfac mystrfac(&myion,&mygrid); @@ -160,6 +164,10 @@ int band_minimizer(MPI_Comm comm_world0, std::string &rtdbstring, std::ostream & /* initialize psps */ CPseudopotential mypsp(&myion,&mygrid,&mystrfac,control,coutput); + // initialize electron operators + cElectron_Operators myelectron(&mygrid,&mykin,&mycoulomb,&myxc,&mypsp); + + /* setup ewald */ Ewald myewald(&myparallel,&myion,&mylattice,control,mypsp.zv); myewald.phafac(); @@ -168,9 +176,9 @@ int band_minimizer(MPI_Comm comm_world0, std::string &rtdbstring, std::ostream & // initialize solid -// Molecule mymolecule(control.input_movecs_filename(), -// control.input_movecs_initialize(),&mygrid,&myion, -// &mystrfac,&myewald,&myelectron,&mypsp,coutput); + Solid mysolid(control.input_movecs_filename(), + control.input_movecs_initialize(),&mygrid,&myion, + &mystrfac,&myewald,&myelectron,&mypsp,coutput); /* intialize the linesearch */ util_linesearch_init(); @@ -388,11 +396,11 @@ int band_minimizer(MPI_Comm comm_world0, std::string &rtdbstring, std::ostream & if (flag < 0) { - //EV = cgsd_noit_energy(mymolecule, true, coutput); + EV = band_cgsd_noit_energy(mysolid, true, coutput); } else { - // EV = cgsd_energy(control, mymolecule, true, coutput); + EV = band_cgsd_energy(control, mysolid, true, coutput); } if (myparallel.is_master()) seconds(&cpu3); diff --git a/Nwpw/nwpwlib/C3dB/CGrid.cpp b/Nwpw/nwpwlib/C3dB/CGrid.cpp index 616272d0..66c82d7b 100644 --- a/Nwpw/nwpwlib/C3dB/CGrid.cpp +++ b/Nwpw/nwpwlib/C3dB/CGrid.cpp @@ -604,7 +604,7 @@ CGrid::CGrid(Parallel *inparall, Lattice *inlattice, int mapping0, int balance0, aqnffts = new (std::nothrow) int[aqmax](); aqnbb = new (std::nothrow) int[aqmax](); atmp = new (std::nothrow) double[2*aqmax*2*nfft3d*nffts_max](); - + bqmax = pfft3_qsize0; if (staged_gpu_fft_pipeline) bqmax += 6; diff --git a/Nwpw/nwpwlib/C3dB/Cneb.cpp b/Nwpw/nwpwlib/C3dB/Cneb.cpp index f158ba01..b1e28b05 100644 --- a/Nwpw/nwpwlib/C3dB/Cneb.cpp +++ b/Nwpw/nwpwlib/C3dB/Cneb.cpp @@ -807,6 +807,32 @@ void Cneb::gh_fftb(double *psi, double *psi_r) } } + +/************************************* + * * + * Cneb::gh_fftb0 * + * * + *************************************/ +void Cneb::gh_fftb0(double *psi, double *psi_r) +{ + + int shift1 = 2*CGrid::npack1_max(); + int shift2 = n2ft3d; + int indx1n = 0; + int indx2n = 0; + for (auto nbq=0; nbq