This thread is locked.Only browsing is available.
Top Page > Browsing
Numerical error of force and total energy by the Gauss-Legendre quadrature in Total_Energy.c
Date: 2012/12/28 13:39
Name: A. M. Ito

Dear Prof. Ozaki

Thank you very much for opening "OpenMX".
I report the bug fix in OpenMX 3.6.1.

In a calculation result, forces acting on atoms have the numerical error which is not small I think. The numerical error is shown as "force(8)" in stdout when "level.of.stdout = 2" in an input file. The numerical error can be checked by judging the forces in all atoms are zero(small) when geometry is stable structure such as BCC and FCC lattices.
The sample of input file to check is as follows:

Atoms.Number 2
Atoms.SpeciesAndCoordinates.Unit Ang # Ang|AU
<Atoms.SpeciesAndCoordinates
1 W 0.7930 0.7930 0.7930 7.0000 7.0000
2 W 2.3790 2.3790 2.3790 7.0000 7.0000
Atoms.SpeciesAndCoordinates>
Atoms.UnitVectors.Unit Ang # Ang|AU
<Atoms.UnitVectors
3.172000 0.000000 0.000000
0.000000 3.172000 0.000000
0.000000 0.000000 3.172000
Atoms.UnitVectors>
scf.EigenvalueSolver band

Here, the cause of this numerical error is the symmetry breaking of integration range of "phi", which means the phi angle in spherical grid, by using Gauss-Legendre quadrature at about 1614-1790 lines in "Total_Energy.c". Just point of the integration of "phi" is 1689-1741 lines. The default integration range of phi is [0, 2.0*PI). However, cos(phi), which is used for x coordinate, is sampled asymmetrically as discrete variables because the Gauss-Legendre quadrature generate unequally-spaced grid in [-1,1] space. In actual, there are much sampling points if the phi is close to 0 and 2.0*PI, while there are few sampling points if the phi is close to PI. This fact causes the numerical difference between the summation of function (functional) cos(phi) and analytical integration.

This problem was not easily solved even if the number of sampling point "Exc0_GL_Mesh2" becomes greater value. However, if is simply solved by the following way. To add symmetric property to sample points of cos(phi), the integration range [0, 2*PI) of phi is divided into [0,PI) and [0,-PI). The Gauss-Legendre quadrature is performed for each range. It is enough that the number of sampling point for each range is half of default Exc0_GL_Mesh2.

This modification easily prevents the numerical error in force(8).
In the next comment, I show my modified code.
メンテ
Page: [1]

Re: Numerical error of force and total energy by the Gauss-Legendre quadrature in Total_Energy.c ( No.1 )
Date: 2012/12/28 13:37
Name: A. M. Ito

The sample code of the above modification is shown here.
The "AITUNE" means my modified point.

--openmx_comment.h, line 51------------------------------------------

//AITUNE//#define Exc0_GL_Mesh2 48 /* # of grids in a Gauss-Legendre quadrature */
#define Exc0_GL_Mesh2 24 /* AITUNE: Gauss-Legendre quadrature range of [0,2*PI) is splitted in [0,PI) and [0,-PI) */


--openmx_comment.h, line 1614------------------------------------------

/****************************************************
calculation of Exc^(0) and its contribution
to forces on the fine mesh
****************************************************/

/* get Nthrds0 */
#pragma omp parallel shared(Nthrds0)
{
Nthrds0 = omp_get_num_threads();
}

/* allocation of arrays */
My_sumr = (double*)malloc(sizeof(double)*Nthrds0);
My_sumrx = (double*)malloc(sizeof(double)*Nthrds0);
My_sumry = (double*)malloc(sizeof(double)*Nthrds0);
My_sumrz = (double*)malloc(sizeof(double)*Nthrds0);

/* start calc. */

rs = 0.0;

ts = 0.0;
te = PI;
St = te + ts;
Dt = te - ts;

ps = 0.0;
//AITUNE// pe = 2.0*PI;
pe = PI;
/* AITUNE: Gauss-Legendre quadrature range of [0,2*PI) is splitted in [0,PI) and [0,-PI). */
/* This modification prevents the numerical error in force(8) which is caused by the */
/* integration for an asymmetric range of "phi" using Gauss-Legendre quadrature. */

Sp = pe + ps;
Dp = pe - ps;

sum = 0.0;

for (Mc_AN=1; Mc_AN<=Matomnum; Mc_AN++){

Gc_AN = M2G[Mc_AN];
Cwan = WhatSpecies[Gc_AN];
re = Spe_Atom_Cut1[Cwan];
Sr = re + rs;
Dr = re - rs;

for (Nloop=0; Nloop<Nthrds0; Nloop++){
My_sumr[Nloop] = 0.0;
My_sumrx[Nloop] = 0.0;
My_sumry[Nloop] = 0.0;
My_sumrz[Nloop] = 0.0;
}

#pragma omp parallel shared(My_sumr,My_sumrx,My_sumry,My_sumrz,Dr,Sr,CoarseGL_Abscissae,CoarseGL_Weight,Dt,St,Exc0_GL_Abscissae1,Dp,Sp,Exc0_GL_Abscissae2,Gxyz,Gc_AN,FNAN,natn,ncn,WhatSpecies,atv,F_Vxc_flag,Cwan,Exc0_GL_Weight2,Exc0_GL_Weight1,PCC_switch) private(OMPID,Nthrds,Nprocs,ir,r,sumt,sumtx,sumty,sumtz,it,theta,sit,cot,sump,sumpx,sumpy,sumpz,ip,phi,x0,y0,z0,h_AN,Gh_AN,Rn,Hwan,x1,y1,z1,dx,dy,dz,r1,den,den0,gden0,dx1,dy1,dz1,exc0,vxc0)
{

/* get info. on OpenMP */

OMPID = omp_get_thread_num();
Nthrds = omp_get_num_threads();
Nprocs = omp_get_num_procs();

for (ir=(OMPID*CoarseGL_Mesh/Nthrds); ir<((OMPID+1)*CoarseGL_Mesh/Nthrds); ir++){

r = 0.50*(Dr*CoarseGL_Abscissae[ir] + Sr);
sumt = 0.0;
sumtx = 0.0;
sumty = 0.0;
sumtz = 0.0;

for (it=0; it<Exc0_GL_Mesh1; it++){

theta = 0.50*(Dt*Exc0_GL_Abscissae1[it] + St);
sit = sin(theta);
cot = cos(theta);
sump = 0.0;
sumpx = 0.0;
sumpy = 0.0;
sumpz = 0.0;

for (ip=0; ip<Exc0_GL_Mesh2; ip++){

phi = 0.50*(Dp*Exc0_GL_Abscissae2[ip] + Sp);
x0 = r*sit*cos(phi) + Gxyz[Gc_AN][1];
y0 = r*sit*sin(phi) + Gxyz[Gc_AN][2];
z0 = r*cot + Gxyz[Gc_AN][3];

/* calculate rho_atom + rho_pcc */

den = 0.0;
for (h_AN=0; h_AN<=FNAN[Gc_AN]; h_AN++){
Gh_AN = natn[Gc_AN][h_AN];
Rn = ncn[Gc_AN][h_AN];
Hwan = WhatSpecies[Gh_AN];

x1 = Gxyz[Gh_AN][1] + atv[Rn][1];
y1 = Gxyz[Gh_AN][2] + atv[Rn][2];
z1 = Gxyz[Gh_AN][3] + atv[Rn][3];

dx = x1 - x0;
dy = y1 - y0;
dz = z1 - z0;
r1 = sqrt(dx*dx + dy*dy + dz*dz);

/* calculate density */

den += ( AtomicDenF(Hwan,r1) + (double)PCC_switch*AtomicPCCF(Hwan,r1) )*F_Vxc_flag;
if (h_AN==0) den0 = den;

/* calculate gradient of density */

if (h_AN==0){
gden0 = ( Dr_AtomicDenF(Cwan,r1) + (double)PCC_switch*Dr_AtomicPCCF(Cwan,r1) )*F_Vxc_flag;

dx1 = gden0/r1*dx;
dy1 = gden0/r1*dy;
dz1 = gden0/r1*dz;
}
}

/* calculate the CA-LDA exchange-correlation energy density */
exc0 = XC_Ceperly_Alder(den,0);

/* calculate the CA-LDA exchange-correlation potential */
vxc0 = XC_Ceperly_Alder(den,1);

/* phi for Gauss-Legendre quadrature */
sump += Exc0_GL_Weight2[ip]*den0*exc0;
sumpx += Exc0_GL_Weight2[ip]*vxc0*dx1;
sumpy += Exc0_GL_Weight2[ip]*vxc0*dy1;
sumpz += Exc0_GL_Weight2[ip]*vxc0*dz1;

} /* ip */

//AITUNE///////////////////////////////////////////////////////////
// The ip-loop is copied to split the integration range of phi. //
// The difference of this second loop from the above one is only //
// that "phi" is changed into "-phi". The others are just copy. //
for (ip=0; ip<Exc0_GL_Mesh2; ip++){

phi = - 0.50*(Dp*Exc0_GL_Abscissae2[ip] + Sp); //AITUNE: different in sign //
x0 = r*sit*cos(phi) + Gxyz[Gc_AN][1];
y0 = r*sit*sin(phi) + Gxyz[Gc_AN][2];
z0 = r*cot + Gxyz[Gc_AN][3];

/* calculate rho_atom + rho_pcc */

den = 0.0;
for (h_AN=0; h_AN<=FNAN[Gc_AN]; h_AN++){
Gh_AN = natn[Gc_AN][h_AN];
Rn = ncn[Gc_AN][h_AN];
Hwan = WhatSpecies[Gh_AN];

x1 = Gxyz[Gh_AN][1] + atv[Rn][1];
y1 = Gxyz[Gh_AN][2] + atv[Rn][2];
z1 = Gxyz[Gh_AN][3] + atv[Rn][3];

dx = x1 - x0;
dy = y1 - y0;
dz = z1 - z0;
r1 = sqrt(dx*dx + dy*dy + dz*dz);

/* calculate density */

den += ( AtomicDenF(Hwan,r1) + (double)PCC_switch*AtomicPCCF(Hwan,r1) )*F_Vxc_flag;
if (h_AN==0) den0 = den;

/* calculate gradient of density */

if (h_AN==0){
gden0 = ( Dr_AtomicDenF(Cwan,r1) + (double)PCC_switch*Dr_AtomicPCCF(Cwan,r1) )*F_Vxc_flag;

dx1 = gden0/r1*dx;
dy1 = gden0/r1*dy;
dz1 = gden0/r1*dz;
}
}

/* calculate the CA-LDA exchange-correlation energy density */
exc0 = XC_Ceperly_Alder(den,0);

/* calculate the CA-LDA exchange-correlation potential */
vxc0 = XC_Ceperly_Alder(den,1);

/* phi for Gauss-Legendre quadrature */
sump += Exc0_GL_Weight2[ip]*den0*exc0;
sumpx += Exc0_GL_Weight2[ip]*vxc0*dx1;
sumpy += Exc0_GL_Weight2[ip]*vxc0*dy1;
sumpz += Exc0_GL_Weight2[ip]*vxc0*dz1;

} /* ip */
////////////////////////////////////AITUNE-fin//

/* theta for Gauss-Legendre quadrature */
sumt += 0.5*Dp*sump*sit*Exc0_GL_Weight1[it];
sumtx += 0.5*Dp*sumpx*sit*Exc0_GL_Weight1[it];
sumty += 0.5*Dp*sumpy*sit*Exc0_GL_Weight1[it];
sumtz += 0.5*Dp*sumpz*sit*Exc0_GL_Weight1[it];

} /* it */

/* r for Gauss-Legendre quadrature */
My_sumr[OMPID] += 0.5*Dt*sumt*r*r*CoarseGL_Weight[ir];
My_sumrx[OMPID] += 0.5*Dt*sumtx*r*r*CoarseGL_Weight[ir];
My_sumry[OMPID] += 0.5*Dt*sumty*r*r*CoarseGL_Weight[ir];
My_sumrz[OMPID] += 0.5*Dt*sumtz*r*r*CoarseGL_Weight[ir];

} /* ir */
} /* #pragma omp */

sumr = 0.0;
for (Nloop=0; Nloop<Nthrds0; Nloop++){
sumr += My_sumr[Nloop];
}

sum += 0.5*Dr*sumr;

/* add force */

sumrx = 0.0;
sumry = 0.0;
sumrz = 0.0;
for (Nloop=0; Nloop<Nthrds0; Nloop++){
sumrx += My_sumrx[Nloop];
sumry += My_sumry[Nloop];
sumrz += My_sumrz[Nloop];
}

if (2<=level_stdout){
printf("<Total_Ene> force(8) myid=%2d Mc_AN=%2d Gc_AN=%2d %15.12f %15.12f %15.12f\n",
myid,Mc_AN,Gc_AN,
0.5*Dr*sumrx, 0.5*Dr*sumry, 0.5*Dr*sumrz);fflush(stdout);
}

Gxyz[Gc_AN][17] += 0.5*Dr*sumrx;
Gxyz[Gc_AN][18] += 0.5*Dr*sumry;
Gxyz[Gc_AN][19] += 0.5*Dr*sumrz;

} /* Mc_AN */

/* add Exc^0 calculated on the fine mesh to My_EXC */

メンテ
Re: Numerical error of force and total energy by the Gauss-Legendre quadrature in Total_Energy.c ( No.2 )
Date: 2012/12/28 13:41
Name: A. M. Ito

Thank you very much for your reading.
Please have a happy new year.

Best reagrds,
Atsushi M. Ito
メンテ

Page: [1]