Skip to content

Commit

Permalink
Merge pull request #188 from StevenZhangCSFM/add_rSCAN_r2SCAN
Browse files Browse the repository at this point in the history
  • Loading branch information
phanish-suryanarayana authored Jul 10, 2023
2 parents 3441f91 + 66ec9de commit 50704ff
Show file tree
Hide file tree
Showing 43 changed files with 16,992 additions and 22 deletions.
6 changes: 6 additions & 0 deletions ChangeLog
Original file line number Diff line number Diff line change
Expand Up @@ -3,6 +3,12 @@
-Name
-changes

--------------
Jul 10, 2023
Name: Boqin Zhang
Changes: (exchangeCorrelation.c, xc/mgga/)
1. Add r2scan metaGGA functional

--------------
Jun 26, 2023
Name: Xin Jing, Boqin Zhang
Expand Down
4 changes: 2 additions & 2 deletions doc/.LaTeX/System.tex
Original file line number Diff line number Diff line change
Expand Up @@ -329,13 +329,13 @@
Choice of exchange-correlation functional. Options are \texttt{LDA\_PW} (Perdew-Wang LDA), \texttt{LDA\_PZ} (Purdew-Zunger LDA),
\texttt{GGA\_PBE} (PBE GGA), \texttt{GGA\_RPBE} (revised PBE GGA), and \texttt{GGA\_PBEsol} (PBE GGA revised for solids),
\texttt{PBE0}, \texttt{HF} (Hartree-Fock), \texttt{HSE},
\texttt{vdWDF1} (van der Waals Density Functional developed by Dion et al.), \texttt{vdWDF2} (vdW Density Functional modified by Lee et al), and \texttt{SCAN} (SCAN metaGGA).
\texttt{vdWDF1} (van der Waals Density Functional developed by Dion et al.), \texttt{vdWDF2} (vdW Density Functional modified by Lee et al), and \texttt{SCAN} (SCAN metaGGA), and \texttt{R2SCAN} (r2SCAN metaGGA).
\end{block}

\begin{block}{Remark}
For spin-polarized calculation (\hyperlink{SPIN_TYP}{\texttt{SPIN\_TYP}} = 1), \texttt{LDA\_PZ} is not available.

Before using \texttt{vdWDF1} or \texttt{vdWDF2}, please read the description and remark of \hyperlink{VDWDF_GEN_KERNEL}{\texttt{VDWDF\_GEN\_KERNEL}} to see the requirements.
Currently SCAN and R2SCAN does not support nonlinear core correction pseudopotential.
\end{block}

\end{frame}
Expand Down
Binary file modified doc/Manual.pdf
Binary file not shown.
26 changes: 23 additions & 3 deletions src/exchangeCorrelation.c
Original file line number Diff line number Diff line change
Expand Up @@ -24,6 +24,7 @@
#include "vdWDFnonlinearCorre.h"
#include "mGGAtauTransferTauVxc.h"
#include "mGGAscan.h"
#include "mGGAr2scan.h"

/**
* @brief Calculate exchange correlation potential
Expand Down Expand Up @@ -98,6 +99,9 @@ void Calculate_Vxc(SPARC_OBJ *pSPARC)
case 4:
scanx(DMnd, rho, sigma, tau, ex, vx, v2x, v3x);
break;
case 6:
r2scanx(DMnd, rho, sigma, tau, ex, vx, v2x, v3x);
break;
default:
memset(ex, 0, sizeof(double) * DMnd);
memset(vx, 0, sizeof(double) * DMnd);
Expand Down Expand Up @@ -127,6 +131,9 @@ void Calculate_Vxc(SPARC_OBJ *pSPARC)
case 4:
scanc(DMnd, rho, sigma, tau, ec, vc, v2c, v3c);
break;
case 6:
r2scanc(DMnd, rho, sigma, tau, ec, vc, v2c, v3c);
break;
default:
memset(ec, 0, sizeof(double) * DMnd);
memset(vc, 0, sizeof(double) * DMnd);
Expand Down Expand Up @@ -248,6 +255,9 @@ void Calculate_Vxc(SPARC_OBJ *pSPARC)
case 4:
scanx_spin(DMnd, rho, sigma, tau, ex, vx, v2x, v3x);
break;
case 6:
r2scanx_spin(DMnd, rho, sigma, tau, ex, vx, v2x, v3x);
break;
default:
memset(ex, 0, sizeof(double) * DMnd);
memset(vx, 0, sizeof(double) * DMnd*2);
Expand Down Expand Up @@ -277,6 +287,9 @@ void Calculate_Vxc(SPARC_OBJ *pSPARC)
case 4:
scanc_spin(DMnd, rho, sigma, tau, ec, vc, v2c, v3c);
break;
case 6:
r2scanc_spin(DMnd, rho, sigma, tau, ec, vc, v2c, v3c);
break;
default:
memset(ec, 0, sizeof(double) * DMnd);
memset(vc, 0, sizeof(double) * DMnd*2);
Expand Down Expand Up @@ -380,9 +393,16 @@ void Calculate_Vxc(SPARC_OBJ *pSPARC)

if (pSPARC->ixc[2]) {
if (pSPARC->countPotentialCalculate == 0) { // restore metaGGA labels after 1st SCF
pSPARC->ixc[0] = 4; pSPARC->ixc[1] = 4;
pSPARC->ixc[2] = 1; pSPARC->ixc[3] = 0;
pSPARC->xcoption[0] = 0; pSPARC->xcoption[1] = 0;
if (strcmpi(pSPARC->XC, "SCAN") == 0) {
pSPARC->ixc[0] = 4; pSPARC->ixc[1] = 4;
pSPARC->ixc[2] = 1; pSPARC->ixc[3] = 0;
pSPARC->xcoption[0] = 0; pSPARC->xcoption[1] = 0;
}
else if (strcmpi(pSPARC->XC, "R2SCAN") == 0) {
pSPARC->ixc[0] = 6; pSPARC->ixc[1] = 6;
pSPARC->ixc[2] = 1; pSPARC->ixc[3] = 0;
pSPARC->xcoption[0] = 0; pSPARC->xcoption[1] = 0;
}
}
}

Expand Down
4 changes: 1 addition & 3 deletions src/finalization.c
Original file line number Diff line number Diff line change
Expand Up @@ -172,9 +172,7 @@ void Free_basic(SPARC_OBJ *pSPARC) {
free(pSPARC->Vc);
free(pSPARC->XCPotential);
free(pSPARC->e_xc);
if(strcmpi(pSPARC->XC,"GGA_PBE") == 0 || strcmpi(pSPARC->XC,"GGA_RPBE") == 0 || strcmpi(pSPARC->XC,"GGA_PBEsol") == 0
|| strcmpi(pSPARC->XC,"PBE0") == 0 || strcmpi(pSPARC->XC,"HF") == 0 || strcmpi(pSPARC->XC,"HSE") == 0 || strcmpi(pSPARC->XC,"SCAN") == 0
|| strcmpi(pSPARC->XC,"vdWDF1") == 0 || strcmpi(pSPARC->XC,"vdWDF2") == 0){
if(pSPARC->isgradient){
free(pSPARC->Dxcdgrho);
}
free(pSPARC->elecstPotential);
Expand Down
11 changes: 8 additions & 3 deletions src/initialization.c
Original file line number Diff line number Diff line change
Expand Up @@ -1247,9 +1247,9 @@ void SPARC_copy_input(SPARC_OBJ *pSPARC, SPARC_INPUT_OBJ *pSPARC_Input) {

// check if exchange-correlation functional is metaGGA
pSPARC->mGGAflag = 0;
if (strcmpi(pSPARC->XC, "SCAN") == 0) { // it can be expand, such as adding r2SCAN
if ((strcmpi(pSPARC->XC, "SCAN") == 0) || (strcmpi(pSPARC->XC, "R2SCAN") == 0)) { // it can be expand, such as adding r2SCAN
if (pSPARC->NLCC_flag) {
if (!rank) printf("\nERROR: currently SCAN functional does not support applying NLCC pseudopotential!\n");
if (!rank) printf("\nERROR: currently SCAN and R2SCAN functional does not support applying NLCC pseudopotential!\n");
exit(EXIT_FAILURE);
}
pSPARC->mGGAflag = 1;
Expand Down Expand Up @@ -3177,7 +3177,7 @@ void write_output_init(SPARC_OBJ *pSPARC) {
}

fprintf(output_fp,"***************************************************************************\n");
fprintf(output_fp,"* SPARC (version Jun 26, 2023) *\n");
fprintf(output_fp,"* SPARC (version Jul 10, 2023) *\n");
fprintf(output_fp,"* Copyright (c) 2020 Material Physics & Mechanics Group, Georgia Tech *\n");
fprintf(output_fp,"* Distributed under GNU General Public License 3 (GPL) *\n");
fprintf(output_fp,"* Start time: %s *\n",c_time_str);
Expand Down Expand Up @@ -4052,6 +4052,11 @@ void xc_decomposition(SPARC_OBJ *pSPARC)
pSPARC->ixc[0] = 4; pSPARC->ixc[1] = 4;
pSPARC->ixc[2] = 1; pSPARC->ixc[3] = 0;
pSPARC->isgradient = 1;
} else if (strcmpi(pSPARC->XC, "R2SCAN") == 0) {
xc = -497498;
pSPARC->ixc[0] = 6; pSPARC->ixc[1] = 6;
pSPARC->ixc[2] = 1; pSPARC->ixc[3] = 0;
pSPARC->isgradient = 1;
} else if (strcmpi(pSPARC->XC, "vdWDF1") == 0) {
xc = -102; // this is the index of Zhang-Yang revPBE exchange in Libxc
pSPARC->ixc[0] = 2; pSPARC->ixc[1] = 2;
Expand Down
2 changes: 1 addition & 1 deletion src/makefile
Original file line number Diff line number Diff line change
Expand Up @@ -89,7 +89,7 @@ OBJSC = main.o initialization.o readfiles.o atomdata.o parallelization.o relax.o
xc/vdW/vdWDF/vdWDFinitialization.o xc/vdW/vdWDF/vdWDFfinalization.o \
xc/vdW/vdWDF/vdWDFnonlinearCorre.o xc/vdW/vdWDF/vdWDFreadKernel.o \
xc/vdW/vdWDF/vdWDFstress.o xc/vdW/vdWDF/vdWDFparallelization.o \
xc/mgga/mGGAscan.o xc/mgga/mGGAhamiltonianTerm.o xc/mgga/mGGAstress.o \
xc/mgga/mGGAscan.o xc/mgga/mGGAr2scan.o xc/mgga/mGGAhamiltonianTerm.o xc/mgga/mGGAstress.o \
xc/mgga/mGGAinitialization.o xc/mgga/mGGAfinalization.o xc/mgga/mGGAtauTransferTauVxc.o \
xc/exx/exactExchange.o xc/exx/exactExchangeKpt.o xc/exx/exactExchangeInitialization.o \
xc/exx/exactExchangeFinalization.o xc/exx/exactExchangeStress.o \
Expand Down
4 changes: 1 addition & 3 deletions src/parallelization.c
Original file line number Diff line number Diff line change
Expand Up @@ -1089,9 +1089,7 @@ void Setup_Comms(SPARC_OBJ *pSPARC) {
assert(pSPARC->e_xc != NULL);

// if GGA then allocate for xc energy per particle for each grid point and der. wrt. grho
if(strcmpi(pSPARC->XC,"GGA_PBE") == 0 || strcmpi(pSPARC->XC,"GGA_RPBE") == 0 || strcmpi(pSPARC->XC,"GGA_PBEsol") == 0
|| strcmpi(pSPARC->XC,"PBE0") == 0 || strcmpi(pSPARC->XC,"HF") == 0 || strcmpi(pSPARC->XC,"HSE") == 0 || strcmpi(pSPARC->XC,"SCAN") == 0
|| strcmpi(pSPARC->XC,"vdWDF1") == 0 || strcmpi(pSPARC->XC,"vdWDF2") == 0){
if(pSPARC->isgradient) {
pSPARC->Dxcdgrho = (double *)malloc( DMnd * (2*pSPARC->Nspin - 1) * sizeof(double) );
assert(pSPARC->Dxcdgrho != NULL);
}
Expand Down
4 changes: 1 addition & 3 deletions src/pressure.c
Original file line number Diff line number Diff line change
Expand Up @@ -142,9 +142,7 @@ void Calculate_XC_pressure(SPARC_OBJ *pSPARC) {

if(strcmpi(pSPARC->XC,"LDA_PW") == 0 || strcmpi(pSPARC->XC,"LDA_PZ") == 0){
pSPARC->pres_xc = 3 * pSPARC->Exc - pSPARC->Exc_corr;
} else if(strcmpi(pSPARC->XC,"GGA_PBE") == 0 || strcmpi(pSPARC->XC,"GGA_RPBE") == 0 || strcmpi(pSPARC->XC,"GGA_PBEsol") == 0
|| strcmpi(pSPARC->XC,"PBE0") == 0 || strcmpi(pSPARC->XC,"HF") == 0 || strcmpi(pSPARC->XC,"HSE") == 0
|| strcmpi(pSPARC->XC,"vdWDF1") == 0 || strcmpi(pSPARC->XC,"vdWDF2") == 0 || strcmpi(pSPARC->XC,"SCAN") == 0){
} else if(pSPARC->isgradient){
pSPARC->pres_xc = 3 * pSPARC->Exc - pSPARC->Exc_corr;
int DMnd, i;
DMnd = (2*pSPARC->Nspin - 1) * pSPARC->Nd_d;
Expand Down
4 changes: 1 addition & 3 deletions src/stress.c
Original file line number Diff line number Diff line change
Expand Up @@ -500,9 +500,7 @@ void Calculate_XC_stress(SPARC_OBJ *pSPARC) {
if(strcmpi(pSPARC->XC,"LDA_PW") == 0 || strcmpi(pSPARC->XC,"LDA_PZ") == 0){
pSPARC->stress_xc[0] = pSPARC->stress_xc[3] = pSPARC->stress_xc[5] = pSPARC->Exc - pSPARC->Exc_corr;
pSPARC->stress_xc[1] = pSPARC->stress_xc[2] = pSPARC->stress_xc[4] = 0.0;
} else if(strcmpi(pSPARC->XC,"GGA_PBE") == 0 || strcmpi(pSPARC->XC,"GGA_RPBE") == 0 || strcmpi(pSPARC->XC,"GGA_PBEsol") == 0
|| strcmpi(pSPARC->XC,"PBE0") == 0 || strcmpi(pSPARC->XC,"HF") == 0 || strcmpi(pSPARC->XC,"HSE") == 0
|| strcmpi(pSPARC->XC,"vdWDF1") == 0 || strcmpi(pSPARC->XC,"vdWDF2") == 0 || strcmpi(pSPARC->XC, "SCAN") == 0){
} else if(pSPARC->isgradient){
pSPARC->stress_xc[0] = pSPARC->stress_xc[3] = pSPARC->stress_xc[5] = pSPARC->Exc - pSPARC->Exc_corr;
pSPARC->stress_xc[1] = pSPARC->stress_xc[2] = pSPARC->stress_xc[4] = 0.0;
int DMnd, i;
Expand Down
162 changes: 162 additions & 0 deletions src/xc/mgga/include/mGGAr2scan.h
Original file line number Diff line number Diff line change
@@ -0,0 +1,162 @@
/**
* @file mGGAr2scan.h
* @brief This file contains the function declarations for SCAN functional.
*
* @authors Boqin Zhang <[email protected]>
* Phanish Suryanarayana <[email protected]>
*
* Copyright (c) 2020 Material Physics & Mechanics Group, Georgia Tech.
*/

#ifndef M_PI
#define M_PI 3.14159265358979323846
#endif

#ifndef R2SCAN_H
#define R2SCAN_H

#include "isddft.h"

/**
* @brief One of the four main functions in mGGAscan.c, compute \epsilon_x and potentials of exchange of SCAN. Called by Calculate_Vxc in exchangeCorrelation.c
*
* @param DMnd number of grid points in the processor in phi domain
* @param rho electron density vector
* @param sigma square of length of gradient of electron density |grad n|^2
* @param tau the kinetic density tau vector
* @param ex the exchange energy density epsilon_x
* @param vx d(n epsilon)/d(n)
* @param v2x d(n epsilon)/d(|grad n|)/|grad n|
* @param v3x d(n epsilon)/d(tau)
*/
void r2scanx(int DMnd, double *rho, double *sigma, double *tau, double *ex, double *vx, double *v2x, double *v3x);

/**
* @brief compute important medium variables: s and ds/dn, ds/d|grad n|; \alpha and d\alpha/dn, d\alpha/d(grad n) and d\alpha/d\tau
*
* @param length length of the rho vector (also the length of |grad n| and tau)
* @param rho electron density vector
* @param normDrho size of gradient of electron density |grad n|
* @param tau the kinetic density tau vector
* @param p_dpdn_dpddn the medium variable p, and dp/dn, dp/d|grad n|
* @param alpha_dadn_daddn_dadtau the medium variable alpha, and d alpha/dn, d alpha/d|grad n|, d alpha/d tau
*/
void basic_r2scanx_variables(int length, double *rho, double *normDrho, double *tau, double **p_dpdn_dpddn, double **alpha_dadn_daddn_dadtau);

/**
* @brief compute exchange part of \epsilon and potential of SCAN
*
* @param length length of the rho vector (also the length of |grad n| and tau)
* @param rho electron density vector
* @param p_dpdn_dpddn the medium variable p, and dp/dn, dp/d|grad n|
* @param alpha_dadn_daddn_dadtau the medium variable alpha, and d alpha/dn, d alpha/d|grad n|, d alpha/d tau
* @param epsilonx exchange energy density
* @param vx1 d(n epsilon_x)/dn
* @param vx2 d(n epsilon_x)/d|grad n|
* @param vx3 d(n epsilon_x)/d tau
*/
void Calculate_r2scanx(int length, double *rho, double **p_dpdn_dpddn, double **alpha_dadn_daddn_dadtau, double *epsilonx, double *vx, double *v2x, double *v3x);

/**
* @brief One of the four main functions in mGGAscan.c, compute \epsilon_c and potentials of correlation of SCAN. Called by Calculate_Vxc in exchangeCorrelation.c
*
* @param DMnd number of grid points in the processor in phi domain
* @param rho electron density vector
* @param sigma square of length of gradient of electron density |grad n|^2
* @param tau the kinetic density tau vector
* @param ec the correlation energy density epsilon_c
* @param vc d(n epsilon)/d(n)
* @param v2c d(n epsilon)/d(|grad n|)/|grad n|
* @param v3c d(n epsilon)/d(tau)
*/
void r2scanc(int DMnd, double *rho, double *sigma, double *tau, double *ec, double *vc, double *v2c, double *v3c);

/**
* @brief compute important medium variables: s and ds/dn, ds/d|grad n|; \alpha and d\alpha/dn, d\alpha/d(grad n) and d\alpha/d\tau
*
* @param length length of the rho vector (also the length of |grad n| and tau)
* @param rho electron density vector
* @param normDrho size of gradient of electron density |grad n|
* @param tau the kinetic density tau vector
* @param s_dsdn_dsddn the medium variable s, and ds/dn, ds/d|grad n|
* @param p_dpdn_dpddn the medium variable p, and dp/dn, dp/d|grad n|
* @param alpha_dadn_daddn_dadtau the medium variable alpha, and d alpha/dn, d alpha/d|grad n|, d alpha/d tau
*/
void basic_r2scanc_variables(int length, double *rho, double *normDrho, double *tau, double **s_dsdn_dsddn, double **p_dpdn_dpddn, double **alpha_dadn_daddn_dadtau);

/**
* @brief compute correlation part of \epsilon and potential of SCAN
*
* @param length length of the rho vector (also the length of |grad n| and tau)
* @param rho electron density vector
* @param s_dsdn_dsddn the medium variable s, and ds/dn, ds/d|grad n|
* @param p_dpdn_dpddn the medium variable p, and dp/dn, dp/d|grad n|
* @param alpha_dadn_daddn_dadtau the medium variable alpha, and d alpha/dn, d alpha/d|grad n|, d alpha/d tau
* @param epsilonc correlation energy density
* @param vc d(n epsilon_c)/dn
* @param v2c d(n epsilon_c)/d|grad n|
* @param v3c d(n epsilon_c)/d tau
*/
void Calculate_r2scanc(int length, double *rho, double **s_dsdn_dsddn, double **p_dpdn_dpddn, double **alpha_dadn_daddn_dadtau, double *epsilonc, double *vc, double *v2c, double *v3c);

/**
* @brief One of the four main functions in mGGAscan.c, compute \epsilon_x and potentials of exchange of SCAN. Called by Calculate_Vxc in exchangeCorrelation.c
*
* @param DMnd number of grid points in the processor in phi domain
* @param rho electron density vector
* @param sigma square of length of gradient of electron density |grad n|^2
* @param tau the kinetic density tau vector
* @param ex the exchange energy density epsilon_x
* @param vx d(n epsilon)/d(n); two vectors: up/dn
* @param v2x d(n epsilon)/d(|grad n|)/|grad n|; two vectors: up/dn
* @param v3x d(n epsilon)/d(tau); two vectors: up/dn
*/
void r2scanx_spin(int DMnd, double *rho, double *sigma, double *tau, double *ex, double *vx, double *v2x, double *v3x);

/**
* @brief One of the four main functions in mGGAscan.c, compute \epsilon_c and potentials of correlation of SCAN. Called by Calculate_Vxc in exchangeCorrelation.c
*
* @param DMnd number of grid points in the processor in phi domain
* @param rho electron density vector
* @param sigma square of length of gradient of electron density |grad n|^2
* @param tau the kinetic density tau vector
* @param ec the correlation energy density epsilon_c
* @param vc d(n epsilon)/d(n), two vectors: up/dn
* @param v2c d(n epsilon)/d(|grad n|)/|grad n|, one vector
* @param v3c d(n epsilon)/d(tau), one vector
*/
void r2scanc_spin(int DMnd, double *rho, double *sigma, double *tau, double *ec, double *vc, double *v2c, double *v3c);

/**
* @brief compute important medium variables: s and ds/dn, ds/d|grad n|; \alpha and d\alpha/dn, d\alpha/d(grad n) and d\alpha/d\tau
*
* @param length length of the rho vector (also the length of |grad n| and tau)
* @param rho electron density vector, n = n_up + n_dn
* @param normDrho length of gradient of electron density |grad n|
* @param tau kinetic density tau vector, tau = tau_up + tau_dn
*
* @param s_dsdn_dsddn the medium variable s, and ds/dn, ds/d|grad n|, from the whole system
* @param p_dpdn_dpddn the medium variable p, and dp/dn, dp/d|grad n|, from the whole system
* @param alpha_dadnup_dadndn_daddn_dadtau the medium variable alpha, and d alpha/dn_up, d alpha/dn_dn, d alpha/d|grad n|, d alpha/d tau
* @param zeta_dzetadnup_dzetadndn the medium varible zeta, and d zeta/dn_up, d zeta/dn_dn, from the whole system
*/
void basic_r2scanc_spin_variables(int length, double *rho, double *normDrho, double *tau, double **s_dsdn_dsddn, double **p_dpdn_dpddn, double **alpha_dadnup_dadndn_daddn_dadtau, double **zeta_dzetadnup_dzetadndn);

/**
* @brief compute correlation part of \epsilon and potential of SCAN
*
* @param length length of the rho vector (also the length of |grad n| and tau)
* @param rho electron density vector, n = n_up + n_dn
* @param s_dsdn_dsddn the medium variable s, and ds/dn, ds/d|grad n|, from the whole system
* @param p_dpdn_dpddn the medium variable p, and dp/dn, dp/d|grad n|, from the whole system
* @param alpha_dadnup_dadndn_daddn_dadtau the medium variable alpha, and d alpha/dn_up, d alpha/dn_dn, d alpha/d|grad n|, d alpha/d tau
* @param zeta_dzetadnup_dzetadndn the medium varible zeta, and d zeta/dn_up, d zeta/dn_dn, from the whole system
*
* @param epsilonc correlation energy density
* @param vc 0~Nd_d-1: d(n epsilon_c)/dn_up; Nd_d~2*Nd_d-1: d(n epsilon_c)/dn_dn
* @param v2c d(n epsilon_c)/d|grad n|
* @param v3c d(n epsilon_c)/d tau
*/
void Calculate_r2scanc_spin(int length, double *rho, double **s_dsdn_dsddn, double **p_dpdn_dpddn, double **alpha_dadnup_dadndn_daddn_dadtau, double **zeta_dzetadnup_dzetadndn, double *epsilonc, double *vc, double *v2c, double *v3c);

#endif // R2SCAN_H
Loading

0 comments on commit 50704ff

Please sign in to comment.