Logo Search packages:      
Sourcecode: openmx version File versions  Download package

Total_Energy.c

/**********************************************************************
  Total_Energy.c:

     Total_Energy.c is a subrutine to calculate the total energy

  Log of Total_Energy.c:

     22/Nov/2001  Released by T.Ozaki
     19/Feb/2006  The subroutine name 'Correction_Energy' was changed 
                  to 'Total_Energy'

***********************************************************************/

#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <time.h>
#include "openmx_common.h"

#ifdef nompi
#include "mimic_mpi.h"
#else
#include "mpi.h"
#endif

static double Calc_Ecore();
static double Calc_EH0(int MD_iter);
static double Calc_Ekin();
static double Calc_Ena();
static double Calc_Enl();
static void Calc_EXC_EH1(double ECE[]);
static double Calc_Ehub();  /* --- added by MJ  */
static void EH0_TwoCenter(int Gc_AN, int h_AN, double VH0ij[4]);
static void EH0_TwoCenter_at_Cutoff(int wan1, int wan2, double VH0ij[4]);


double Total_Energy(int MD_iter, double *****CDM, double ECE[])
{ 
  double time0;
  double TStime,TEtime;
  int numprocs,myid;

  /* MPI */
  MPI_Comm_size(mpi_comm_level1,&numprocs);
  MPI_Comm_rank(mpi_comm_level1,&myid);

  dtime(&TStime);

  /****************************************************
               core-core repulsion energy
  ****************************************************/

  ECE[0] = Calc_Ecore();
  
  /****************************************************
              EH0 = -1/2\int n^a(r) V^a_H dr
  ****************************************************/

  ECE[1] = Calc_EH0(MD_iter);

  /****************************************************
                    kinetic energy
  ****************************************************/

  if (F_Kin_flag==1)  ECE[2] = Calc_Ekin();

  /****************************************************
              neutral atom potential energy
  ****************************************************/

  if (F_VNA_flag==1)  ECE[3] = Calc_Ena();

  /****************************************************
           non-local pseudo potential energy
  ****************************************************/

  if (F_NL_flag==1)  ECE[4] = Calc_Enl();

  /****************************************************
     EXC = \sum_{\sigma}
           n_{\sigma}(\epsilon_{xc}-\mu_{xc,\sigma})

     EH1 = -1/2\int {n(r)+n^a(r)} \delta V_H dr
  ****************************************************/

  Calc_EXC_EH1(ECE);

  /****************************************************
    LDA+U energy   --- added by MJ
  ****************************************************/

  if (F_U_flag==1){
    if (Hub_U_switch==1)  ECE[8] = Calc_Ehub();
    else                  ECE[8] = 0.0;
  }

  /* computational time */

  dtime(&TEtime);
  time0 = TEtime - TStime;
  return time0;
}








double Calc_Ekin()
{
  /****************************************************
          calculate the kinetic energy, Ekin
  ****************************************************/

  int i,j,spin;
  int Mc_AN,Gc_AN,Cwan,Rn,h_AN,Gh_AN,Hwan;
  double My_Ekin,Ekin,Zc,Zh,dum,dum2;
  int numprocs,myid;
  double Stime_atom, Etime_atom;

  /* MPI */
  MPI_Comm_size(mpi_comm_level1,&numprocs);
  MPI_Comm_rank(mpi_comm_level1,&myid);

  My_Ekin = 0.0;

  if (SpinP_switch==0 || SpinP_switch==1){

    for (spin=0; spin<=SpinP_switch; spin++){
      for (Mc_AN=1; Mc_AN<=Matomnum; Mc_AN++){

      dtime(&Stime_atom);

      Gc_AN = M2G[Mc_AN];
      Cwan = WhatSpecies[Gc_AN];

        if (Cnt_switch==0){
        for (h_AN=0; h_AN<=FNAN[Gc_AN]; h_AN++){
          Gh_AN = natn[Gc_AN][h_AN];
          Hwan = WhatSpecies[Gh_AN];
          for (i=0; i<Spe_Total_CNO[Cwan]; i++){
            for (j=0; j<Spe_Total_CNO[Hwan]; j++){
            My_Ekin += DM[0][spin][Mc_AN][h_AN][i][j]*H0[0][Mc_AN][h_AN][i][j];
            }
          }
        }
      }

        else{
        for (h_AN=0; h_AN<=FNAN[Gc_AN]; h_AN++){
          Gh_AN = natn[Gc_AN][h_AN];
          Hwan = WhatSpecies[Gh_AN];
          for (i=0; i<Spe_Total_CNO[Cwan]; i++){
            for (j=0; j<Spe_Total_CNO[Hwan]; j++){
            My_Ekin += DM[0][spin][Mc_AN][h_AN][i][j]*CntH0[0][Mc_AN][h_AN][i][j];
            }
          }
        }
        }

        dtime(&Etime_atom);
        time_per_atom[Gc_AN] += Etime_atom - Stime_atom;
      }
    }

    if (SpinP_switch==0) My_Ekin = 2.0*My_Ekin;

  }
  else if (SpinP_switch==3){

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

      dtime(&Stime_atom);

      Gc_AN = M2G[Mc_AN];
      Cwan = WhatSpecies[Gc_AN];

      for (h_AN=0; h_AN<=FNAN[Gc_AN]; h_AN++){
      Gh_AN = natn[Gc_AN][h_AN];
      Hwan = WhatSpecies[Gh_AN];
      for (i=0; i<Spe_Total_CNO[Cwan]; i++){
        for (j=0; j<Spe_Total_CNO[Hwan]; j++){
          My_Ekin += (DM[0][0][Mc_AN][h_AN][i][j] + DM[0][1][Mc_AN][h_AN][i][j])*H0[0][Mc_AN][h_AN][i][j];
        }
      }
      }

      dtime(&Etime_atom);
      time_per_atom[Gc_AN] += Etime_atom - Stime_atom;
    }
  }

  MPI_Allreduce(&My_Ekin, &Ekin, 1, MPI_DOUBLE, MPI_SUM, mpi_comm_level1);

  return Ekin;  
}





double Calc_Ena()
{
  /****************************************************
     calculate the neutral atom potential energy, Ena
  ****************************************************/

  int i,j,spin;
  int Mc_AN,Gc_AN,Cwan,Rn,h_AN,Gh_AN,Hwan;
  double My_Ena,Ena,Zc,Zh,dum,dum2;
  int numprocs,myid;
  double Stime_atom, Etime_atom;

  /* MPI */
  MPI_Comm_size(mpi_comm_level1,&numprocs);
  MPI_Comm_rank(mpi_comm_level1,&myid);

  My_Ena = 0.0;

  if (SpinP_switch==0 || SpinP_switch==1){

    if (Cnt_switch==1){
      /* temporaly, borrow the CntH0 matrix */
      Cont_Matrix0(HVNA,CntH0[0]);
    }

    for (spin=0; spin<=SpinP_switch; spin++){
      for (Mc_AN=1; Mc_AN<=Matomnum; Mc_AN++){

      dtime(&Stime_atom);

      Gc_AN = M2G[Mc_AN];
      Cwan = WhatSpecies[Gc_AN];

        if (Cnt_switch==0){
        for (h_AN=0; h_AN<=FNAN[Gc_AN]; h_AN++){
          Gh_AN = natn[Gc_AN][h_AN];
          Hwan = WhatSpecies[Gh_AN];
          for (i=0; i<Spe_Total_CNO[Cwan]; i++){
            for (j=0; j<Spe_Total_CNO[Hwan]; j++){
            My_Ena += DM[0][spin][Mc_AN][h_AN][i][j]*HVNA[Mc_AN][h_AN][i][j];
            }
          }
        }
      }

        else {
        for (h_AN=0; h_AN<=FNAN[Gc_AN]; h_AN++){
          Gh_AN = natn[Gc_AN][h_AN];
          Hwan = WhatSpecies[Gh_AN];
          for (i=0; i<Spe_Total_CNO[Cwan]; i++){
            for (j=0; j<Spe_Total_CNO[Hwan]; j++){
            My_Ena += DM[0][spin][Mc_AN][h_AN][i][j]*CntH0[0][Mc_AN][h_AN][i][j];
            }
          }
        }
      }

        dtime(&Etime_atom);
        time_per_atom[Gc_AN] += Etime_atom - Stime_atom;
      }
    }

    if (SpinP_switch==0) My_Ena = 2.0*My_Ena;

  }
  else if (SpinP_switch==3){

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

      dtime(&Stime_atom);

      Gc_AN = M2G[Mc_AN];
      Cwan = WhatSpecies[Gc_AN];

      for (h_AN=0; h_AN<=FNAN[Gc_AN]; h_AN++){
      Gh_AN = natn[Gc_AN][h_AN];
      Hwan = WhatSpecies[Gh_AN];
      for (i=0; i<Spe_Total_CNO[Cwan]; i++){
        for (j=0; j<Spe_Total_CNO[Hwan]; j++){
          My_Ena += (DM[0][0][Mc_AN][h_AN][i][j]+DM[0][1][Mc_AN][h_AN][i][j])*HVNA[Mc_AN][h_AN][i][j];
        }
      }
      }

      dtime(&Etime_atom);
      time_per_atom[Gc_AN] += Etime_atom - Stime_atom;
    }
  }

  MPI_Allreduce(&My_Ena, &Ena, 1, MPI_DOUBLE, MPI_SUM, mpi_comm_level1);

  return Ena;  
}





double Calc_Enl()
{
  /****************************************************
     calculate the non-local pseudo potential energy
  ****************************************************/

  int i,j,spin;
  int Mc_AN,Gc_AN,Cwan,Rn,h_AN,Gh_AN,Hwan;
  double My_Enl,Enl,Zc,Zh,dum,dum2;
  int numprocs,myid;
  double Stime_atom, Etime_atom;

  /* MPI */
  MPI_Comm_size(mpi_comm_level1,&numprocs);
  MPI_Comm_rank(mpi_comm_level1,&myid);

  My_Enl = 0.0;

  if (SpinP_switch==0 || SpinP_switch==1){

    for (spin=0; spin<=SpinP_switch; spin++){

      if (Cnt_switch==1){
        /* temporaly, borrow the CntH0 matrix */
        Cont_Matrix0(HNL[spin],CntH0[0]);
      }

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

      dtime(&Stime_atom);

      Gc_AN = M2G[Mc_AN];
      Cwan = WhatSpecies[Gc_AN];

        if (Cnt_switch==0){

        for (h_AN=0; h_AN<=FNAN[Gc_AN]; h_AN++){
          Gh_AN = natn[Gc_AN][h_AN];
          Hwan = WhatSpecies[Gh_AN];
          for (i=0; i<Spe_Total_CNO[Cwan]; i++){
            for (j=0; j<Spe_Total_CNO[Hwan]; j++){
            My_Enl += DM[0][spin][Mc_AN][h_AN][i][j]*HNL[spin][Mc_AN][h_AN][i][j];
            }
          }
        }
      }

        else {
        for (h_AN=0; h_AN<=FNAN[Gc_AN]; h_AN++){
          Gh_AN = natn[Gc_AN][h_AN];
          Hwan = WhatSpecies[Gh_AN];
          for (i=0; i<Spe_Total_CNO[Cwan]; i++){
            for (j=0; j<Spe_Total_CNO[Hwan]; j++){
            My_Enl += DM[0][spin][Mc_AN][h_AN][i][j]*CntH0[0][Mc_AN][h_AN][i][j];
            }
          }
        }
      }

        dtime(&Etime_atom);
        time_per_atom[Gc_AN] += Etime_atom - Stime_atom;
      }
    }

    if (SpinP_switch==0) My_Enl = 2.0*My_Enl;

  }
  else if (SpinP_switch==3){

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

      dtime(&Stime_atom);

      Gc_AN = M2G[Mc_AN];
      Cwan = WhatSpecies[Gc_AN];

      for (h_AN=0; h_AN<=FNAN[Gc_AN]; h_AN++){
      Gh_AN = natn[Gc_AN][h_AN];
      Hwan = WhatSpecies[Gh_AN];
      for (i=0; i<Spe_Total_CNO[Cwan]; i++){
        for (j=0; j<Spe_Total_CNO[Hwan]; j++){

          My_Enl += 
                 DM[0][0][Mc_AN][h_AN][i][j]*  HNL[0][Mc_AN][h_AN][i][j]
            - iDM[0][0][Mc_AN][h_AN][i][j]*iHNL0[0][Mc_AN][h_AN][i][j]
            +  DM[0][1][Mc_AN][h_AN][i][j]*  HNL[1][Mc_AN][h_AN][i][j]
            - iDM[0][1][Mc_AN][h_AN][i][j]*iHNL0[1][Mc_AN][h_AN][i][j]
         + 2.0*DM[0][2][Mc_AN][h_AN][i][j]*  HNL[2][Mc_AN][h_AN][i][j] 
         - 2.0*DM[0][3][Mc_AN][h_AN][i][j]*iHNL0[2][Mc_AN][h_AN][i][j];
 
        }
      }
      }


      dtime(&Etime_atom);
      time_per_atom[Gc_AN] += Etime_atom - Stime_atom;
    }
  }

  MPI_Allreduce(&My_Enl, &Enl, 1, MPI_DOUBLE, MPI_SUM, mpi_comm_level1);

  return Enl;  
}






double Calc_Ecore()
{
  /****************************************************
                         Ecore
  ****************************************************/

  int Mc_AN,Gc_AN,Cwan,Rn,h_AN,Gh_AN,Hwan;
  double My_Ecore,Ecore,Zc,Zh,dum,dum2;
  double dEx,dEy,dEz,r,lx,ly,lz;
  int numprocs,myid;
  double Stime_atom, Etime_atom;

  /* MPI */
  MPI_Comm_size(mpi_comm_level1,&numprocs);
  MPI_Comm_rank(mpi_comm_level1,&myid);

  if (myid==Host_ID)  printf("  Force calculation #6\n");fflush(stdout);

  My_Ecore = 0.0;

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

    dtime(&Stime_atom);

    Gc_AN = M2G[Mc_AN];
    Cwan = WhatSpecies[Gc_AN];
    Zc = Spe_Core_Charge[Cwan];
    dEx = 0.0;
    dEy = 0.0;
    dEz = 0.0;
    for (h_AN=1; h_AN<=FNAN[Gc_AN]; h_AN++){
      Gh_AN = natn[Gc_AN][h_AN];
      Rn = ncn[Gc_AN][h_AN];
      Hwan = WhatSpecies[Gh_AN];
      Zh = Spe_Core_Charge[Hwan];
      r = Dis[Gc_AN][h_AN];

      /* for empty atoms or finite elemens basis */
      if (r<1.0e-10) r = 1.0e-10;

      lx = (Gxyz[Gc_AN][1] - Gxyz[Gh_AN][1] - atv[Rn][1])/r;
      ly = (Gxyz[Gc_AN][2] - Gxyz[Gh_AN][2] - atv[Rn][2])/r;
      lz = (Gxyz[Gc_AN][3] - Gxyz[Gh_AN][3] - atv[Rn][3])/r;
      dum = Zc*Zh/r;
      dum2 = dum/r;
      My_Ecore = My_Ecore + dum;
      dEx = dEx - lx*dum2;
      dEy = dEy - ly*dum2;
      dEz = dEz - lz*dum2;
    }

    /****************************************************
                        #6 of force
         Contribution from the core-core repulsions
    ****************************************************/

    Gxyz[Gc_AN][17] = Gxyz[Gc_AN][17] + dEx;
    Gxyz[Gc_AN][18] = Gxyz[Gc_AN][18] + dEy;
    Gxyz[Gc_AN][19] = Gxyz[Gc_AN][19] + dEz;

    if (2<=level_stdout){
      printf("<Total_Ene>  force(6) myid=%2d  Mc_AN=%2d Gc_AN=%2d  %15.12f %15.12f %15.12f\n",
             myid,Mc_AN,Gc_AN,dEx,dEy,dEz);fflush(stdout);
    }

    dtime(&Etime_atom);
    time_per_atom[Gc_AN] += Etime_atom - Stime_atom;
  }
  My_Ecore = 0.50*My_Ecore;

  MPI_Allreduce(&My_Ecore, &Ecore, 1, MPI_DOUBLE, MPI_SUM, mpi_comm_level1);

  return Ecore;  
}


double Calc_EH0(int MD_iter)
{
  /****************************************************
              EH0 = -1/2\int n^a(r) V^a_H dr
  ****************************************************/

  int Mc_AN,Gc_AN,h_AN,Gh_AN,num,gnum;
  int wan,wan1,wan2,Nd,n1,n2,n3;
  double bc,dx,x,y,z,r1,r2,rho0;
  double Rcut12,Scale_Grid_Ecut;
  double EH0ij[4],My_EH0,EH0,tmp0;
  double *Fx,*Fy,*Fz,*g0;
  double dEx,dEy,dEz;
  double Z1,Z2,factor;
  double My_dEx,My_dEy,My_dEz;
  int numprocs,myid,ID;
  double Stime_atom, Etime_atom;

  /* MPI */
  MPI_Comm_size(mpi_comm_level1,&numprocs);
  MPI_Comm_rank(mpi_comm_level1,&myid);

  /****************************************************
     allocation of arrays:

    double Fx[Matomnum+1];
    double Fy[Matomnum+1];
    doubel Fz[Matomnum+1];
  ****************************************************/

  Fx = (double*)malloc(sizeof(double)*(Matomnum+1));
  Fy = (double*)malloc(sizeof(double)*(Matomnum+1));
  Fz = (double*)malloc(sizeof(double)*(Matomnum+1));

  /****************************************************
             Set of atomic density on grids
  ****************************************************/

  if (MD_iter==1){

    Scale_Grid_Ecut = 2.0;

    /* estimate sizes of an array g0 */

    Max_Nd = 0;
    num = -1;
    for (wan=0; wan<SpeciesNum; wan++){
      num++;
      Spe_Spe2Ban[wan] = num;
      bc = Spe_Atom_Cut1[wan];
      dx = PI/sqrt(Grid_Ecut*Scale_Grid_Ecut);
      Nd = 2*(int)(bc/dx);
      if (Max_Nd<Nd) Max_Nd = Nd;
    }

    /* allocation of arrays g0 */

    g0 = (double*)malloc(sizeof(double)*Max_Nd);

    /* estimate sizes of arrays GridX,Y,Z_EH0, Arho_EH0, and Wt_EH0 */

    Max_TGN_EH0 = 0;
    num = -1;
    for (wan=0; wan<SpeciesNum; wan++){
      num++;
      Spe_Spe2Ban[wan] = num;
      bc = Spe_Atom_Cut1[wan];
      dx = PI/sqrt(Grid_Ecut*Scale_Grid_Ecut);
      Nd = 2*(int)(bc/dx);
      dx = 2.0*bc/(double)Nd;

      for (n1=0; n1<Nd; n1++){
        g0[n1] = dx*(double)n1 - bc;
      }

      Rcut12 = Spe_Atom_Cut1[wan]*Spe_Atom_Cut1[wan];
      
      gnum = -1; 
      for (n1=Nd/2; n1<Nd; n1++){
      x = g0[n1];
      for (n2=Nd/2; n2<=n1; n2++){
        y = g0[n2];
        tmp0 = x*x + y*y;
        for (n3=0; n3<Nd; n3++){
          z = g0[n3];
          r1 = tmp0 + z*z;
          if (r1<=Rcut12){
            gnum++;  
          }
        }
      }
      }
      
      if (Max_TGN_EH0<(gnum+1)) Max_TGN_EH0 = gnum + 1;

      if (2<=level_stdout){
        printf("<Calc_EH0> A spe=%2d 1D-grids=%2d 3D-grids=%2d\n",wan,Nd,gnum+1);fflush(stdout);
      }
    }

    /* allocation of arrays GridX,Y,Z_EH0, Arho_EH0, and Wt_EH0 */

    Max_TGN_EH0 += 10; 
    Allocate_Arrays(4);

    /* calculate GridX,Y,Z_EH0 and Wt_EH0 */

    num = -1;
    for (wan=0; wan<SpeciesNum; wan++){

      num++;
      Spe_Spe2Ban[wan] = num;
      bc = Spe_Atom_Cut1[wan];
      dx = PI/sqrt(Grid_Ecut*Scale_Grid_Ecut);
      Nd = 2*(int)(bc/dx);
      dx = 2.0*bc/(double)Nd;
      dv_EH0[num] = dx*dx*dx;

      for (n1=0; n1<Nd; n1++){
        g0[n1] = dx*(double)n1 - bc;
      }

      Rcut12 = Spe_Atom_Cut1[wan]*Spe_Atom_Cut1[wan];

      gnum = -1; 
      for (n1=Nd/2; n1<Nd; n1++){
      x = g0[n1];
      for (n2=Nd/2; n2<=n1; n2++){
        y = g0[n2];
        tmp0 = x*x + y*y;
        for (n3=0; n3<Nd; n3++){
          z = g0[n3];
          r1 = tmp0 + z*z;
          if (r1<=Rcut12){

            gnum++;  

              r1 = sqrt(r1);

            GridX_EH0[num][gnum] = x;
            GridY_EH0[num][gnum] = y;
            GridZ_EH0[num][gnum] = z;

            rho0 = AtomicDenF(wan,r1);
            Arho_EH0[num][gnum] = rho0;

            if (n1==Nd/2 && n2==Nd/2)
            Wt_EH0[num][gnum] = 1.0;
            else if (n2==Nd/2 || n1==n2)
            Wt_EH0[num][gnum] = 4.0;
            else 
            Wt_EH0[num][gnum] = 8.0;

          }
        }
      }
      }

      TGN_EH0[num] = gnum + 1;

      if (2<=level_stdout){
        printf("<Calc_EH0> B spe=%2d 1D-grids=%2d 3D-grids=%2d\n",wan,Nd,gnum+1);fflush(stdout);
      }

    }

    /* free */
    free(g0);
  }

  /****************************************************
    calculation of scaling factors:
  ****************************************************/

  if (MD_iter==1){

    for (wan1=0; wan1<SpeciesNum; wan1++){
      r1 = Spe_Atom_Cut1[wan1];
      Z1 = Spe_Core_Charge[wan1];
      for (wan2=0; wan2<SpeciesNum; wan2++){
        EH0_TwoCenter_at_Cutoff(wan1, wan2, EH0ij);
        r2 = Spe_Atom_Cut1[wan2];
        Z2 = Spe_Core_Charge[wan2];
        tmp0 = Z1*Z2/(r1+r2);

        if (1.0e-20<fabs(EH0ij[0])){ 
          EH0_scaling[wan1][wan2] = tmp0/EH0ij[0];
      }
        else{
          EH0_scaling[wan1][wan2] = 0.0;
        }

      }
    }
  }

  /****************************************************
                -1/2\int n^a(r) V^a_H dr
  ****************************************************/

  for (Mc_AN=1; Mc_AN<=Matomnum; Mc_AN++){
    Fx[Mc_AN] = 0.0;
    Fy[Mc_AN] = 0.0;
    Fz[Mc_AN] = 0.0;
  }

  My_EH0 = 0.0;

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

    dtime(&Stime_atom);

    Gc_AN = M2G[Mc_AN];
    wan1 = WhatSpecies[Gc_AN];

    for (h_AN=0; h_AN<=FNAN[Gc_AN]; h_AN++){
      Gh_AN = natn[Gc_AN][h_AN];
      wan2 = WhatSpecies[Gh_AN];

      if (h_AN==0) factor = 1.0;
      else         factor = EH0_scaling[wan1][wan2];
      
      EH0_TwoCenter(Gc_AN, h_AN, EH0ij);
      My_EH0 -= 0.250*factor*EH0ij[0];
      Fx[Mc_AN] = Fx[Mc_AN] - 0.5*factor*EH0ij[1];
      Fy[Mc_AN] = Fy[Mc_AN] - 0.5*factor*EH0ij[2];
      Fz[Mc_AN] = Fz[Mc_AN] - 0.5*factor*EH0ij[3];

      if (h_AN==0) factor = 1.0;
      else         factor = EH0_scaling[wan2][wan1];

      EH0_TwoCenter(Gh_AN, RMI1[Mc_AN][h_AN][0], EH0ij);
      My_EH0 -= 0.250*factor*EH0ij[0];
      Fx[Mc_AN] = Fx[Mc_AN] + 0.5*factor*EH0ij[1];
      Fy[Mc_AN] = Fy[Mc_AN] + 0.5*factor*EH0ij[2];
      Fz[Mc_AN] = Fz[Mc_AN] + 0.5*factor*EH0ij[3];

    }

    dtime(&Etime_atom);
    time_per_atom[Gc_AN] += Etime_atom - Stime_atom;
  }

  MPI_Allreduce(&My_EH0, &EH0, 1, MPI_DOUBLE, MPI_SUM, mpi_comm_level1);

  /****************************************************
                      #7 of force
     Contribution from the Hartree energy between
   the neutral atomic charge and the neutral potential 
  ****************************************************/

  if (myid==Host_ID)  printf("  Force calculation #7\n");fflush(stdout);

  for (Mc_AN=1; Mc_AN<=Matomnum; Mc_AN++){
    Gc_AN = M2G[Mc_AN];

    if (2<=level_stdout){
      printf("<Total_Ene>  force(7) myid=%2d  Mc_AN=%2d Gc_AN=%2d  %15.12f %15.12f %15.12f\n",
              myid,Mc_AN,Gc_AN,Fx[Mc_AN],Fy[Mc_AN],Fz[Mc_AN]);fflush(stdout);
    }

    Gxyz[Gc_AN][17] += Fx[Mc_AN];
    Gxyz[Gc_AN][18] += Fy[Mc_AN];
    Gxyz[Gc_AN][19] += Fz[Mc_AN];
  }

  /****************************************************
                   calculated forces
  ****************************************************/

  My_dEx = 0.0;
  My_dEy = 0.0;
  My_dEz = 0.0;

  for (Mc_AN=1; Mc_AN<=Matomnum; Mc_AN++){
    Gc_AN = M2G[Mc_AN];

    if (2<=level_stdout){
      printf("<Total_Ene>  force(t) myid=%2d  Mc_AN=%2d Gc_AN=%2d  %15.12f %15.12f %15.12f\n",
              myid,Mc_AN,Gc_AN,Gxyz[Gc_AN][17],Gxyz[Gc_AN][18],Gxyz[Gc_AN][19]);fflush(stdout);
    }

    My_dEx += Gxyz[Gc_AN][17];
    My_dEy += Gxyz[Gc_AN][18];
    My_dEz += Gxyz[Gc_AN][19];
  }
  
  MPI_Allreduce(&My_dEx, &dEx, 1, MPI_DOUBLE, MPI_SUM, mpi_comm_level1);
  MPI_Allreduce(&My_dEy, &dEy, 1, MPI_DOUBLE, MPI_SUM, mpi_comm_level1);
  MPI_Allreduce(&My_dEz, &dEz, 1, MPI_DOUBLE, MPI_SUM, mpi_comm_level1);

  /*
  printf("dEx = %18.15f\n",dEx);
  printf("dEy = %18.15f\n",dEy);
  printf("dEz = %18.15f\n",dEz);

  MPI_Finalize();
  exit(1);
  */

  dEx = dEx/(double)atomnum;
  dEy = dEy/(double)atomnum;
  dEz = dEz/(double)atomnum;

  /****************************************************
     correct forces so that the action and reaction
                  raw is satisfied
  ****************************************************/

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

    Gc_AN = M2G[Mc_AN];
    Gxyz[Gc_AN][17] = Gxyz[Gc_AN][17] - dEx;
    Gxyz[Gc_AN][18] = Gxyz[Gc_AN][18] - dEy;
    Gxyz[Gc_AN][19] = Gxyz[Gc_AN][19] - dEz;

    if (2<=level_stdout){
      printf("<Total_Ene>  force(c) myid=%2d  Mc_AN=%2d Gc_AN=%2d  %15.12f %15.12f %15.12f\n",
              myid,Mc_AN,Gc_AN,Gxyz[Gc_AN][17],Gxyz[Gc_AN][18],Gxyz[Gc_AN][19]);
    }
  }
  */

  /****************************************************
   MPI, Gxyz[Gc_AN][17-19]
  ****************************************************/

  for (Gc_AN=1; Gc_AN<=atomnum; Gc_AN++){
    ID = G2ID[Gc_AN];
    MPI_Bcast(&Gxyz[Gc_AN][17], 3, MPI_DOUBLE, ID, mpi_comm_level1);
  }

  /****************************************************
    freeing of arrays:

    double Fx[Matomnum+1];
    double Fy[Matomnum+1];
    doubel Fz[Matomnum+1];
  ****************************************************/

  free(Fx);
  free(Fy);
  free(Fz);

  /* return */

  return EH0;  
}


void Calc_EXC_EH1(double ECE[])
{
  /****************************************************
     EXC = \sum_{\sigma}
           n_{\sigma}(\epsilon_{xc}-\mu_{xc,\sigma})

     EH1 = -1/2\int {n(r)+n^a(r)} \delta V_H dr
  ****************************************************/

  int i,spin,XC_P_switch,max_spin;
  int numS,numR,My_GNum;
  int MN,MN0,MN1,MN2,n1,n2,n3;
  int MNN0,MNN1,nn0,nn1;
  double EXC[2],EH1,sum;
  double My_EXC[2],My_EH1;
  double My_Eef,Eef;
  double sum_charge,My_charge;
  double *tmp_array0,*tmp_array1;
  double **Den,*ADen,*PDen;
  double *dVHart,**Vxc;
  double *Vef;
  int numprocs,myid,tag=999,ID,IDS,IDR;
  double Stime_atom, Etime_atom;

  /* dipole moment */
  int Gc_AN,Mc_AN,spe;
  double x,y,z,den,charge,cden_BG;
  double E_dpx,E_dpy,E_dpz; 
  double E_dpx_BG,E_dpy_BG,E_dpz_BG; 
  double C_dpx,C_dpy,C_dpz;
  double My_E_dpx_BG,My_E_dpy_BG,My_E_dpz_BG; 
  double My_E_dpx,My_E_dpy,My_E_dpz; 
  double My_C_dpx,My_C_dpy,My_C_dpz;
  double AU2Debye,AbsD;
  char file_DPM[YOUSO10] = ".dpm";
  FILE *fp_DPM;
  char buf[fp_bsize];          /* setvbuf */
  MPI_Status stat;
  MPI_Request request;

  /* MPI */
  MPI_Comm_size(mpi_comm_level1,&numprocs);
  MPI_Comm_rank(mpi_comm_level1,&myid);

  /****************************************************
   set Vxc_Grid
  ****************************************************/

  XC_P_switch = 0;
  Set_XC_Grid(XC_P_switch,XC_switch);

  /****************************************************
   MPI:

   Density_Grid[0]
   Density_Grid[1]
   ADensity_Grid
   PCCDensity_Grid
   dVHart_Grid
   Vxc_Grid[0]
   Vxc_Grid[1]
   Vef
  ****************************************************/
  
  My_GNum = My_NGrid1_Poisson*Ngrid2*Ngrid3;  

  Den = (double**)malloc(sizeof(double*)*2); 
  Den[0] = (double*)malloc(sizeof(double)*My_GNum); 
  Den[1] = (double*)malloc(sizeof(double)*My_GNum); 
  ADen = (double*)malloc(sizeof(double)*My_GNum); 
  PDen = (double*)malloc(sizeof(double)*My_GNum); 
  dVHart = (double*)malloc(sizeof(double)*My_GNum); 
  Vxc    = (double**)malloc(sizeof(double*)*2);
  Vxc[0] = (double*)malloc(sizeof(double)*My_GNum); 
  Vxc[1] = (double*)malloc(sizeof(double)*My_GNum);
  Vef = (double*)malloc(sizeof(double)*My_GNum);

  /* initialize */
  for (MN1=0; MN1<My_GNum; MN1++){
    Den[0][MN1] = 0.0;
    Den[1][MN1] = 0.0;
    ADen[MN1]   = 0.0;
    PDen[MN1]   = 0.0;
    dVHart[MN1] = 0.0;
    Vxc[0][MN1] = 0.0;
    Vxc[1][MN1] = 0.0;
    Vef[MN1]    = 0.0;
  }

  /* use their densities using MPI  */
  for (ID=0; ID<numprocs; ID++){

    IDS = (myid + ID) % numprocs;
    IDR = (myid - ID + numprocs) % numprocs;

    /* Isend */
    if (Num_Snd_Grid1[IDS]!=0){

      numS = Num_Snd_Grid1[IDS]*Ngrid2*Ngrid3;        
      tmp_array0 = (double*)malloc(sizeof(double)*numS*8);

      for (i=0; i<Num_Snd_Grid1[IDS]; i++){ 
      n1 = Snd_Grid1[IDS][i];
        nn1 = My_Cell0[n1];

        MN1 = nn1*Ngrid2*Ngrid3;
        MN0 = i*Ngrid2*Ngrid3;
        for (n2=0; n2<Ngrid2; n2++){
          MN2 = n2*Ngrid3;
          for (n3=0; n3<Ngrid3; n3++){
            MN = MN1 + MN2 + n3;
            MNN0 = MN0 + MN2 + n3; 
            tmp_array0[         MNN0] = Density_Grid[0][MN];
            tmp_array0[  numS + MNN0] = Density_Grid[1][MN];
            tmp_array0[2*numS + MNN0] = ADensity_Grid[MN];
            tmp_array0[3*numS + MNN0] = PCCDensity_Grid[MN];
            tmp_array0[4*numS + MNN0] = dVHart_Grid[MN];
            tmp_array0[5*numS + MNN0] = Vxc_Grid[0][MN];
            tmp_array0[6*numS + MNN0] = Vxc_Grid[1][MN];
            tmp_array0[7*numS + MNN0] = VEF_Grid[MN];
        }
      }
      }

      MPI_Isend(&tmp_array0[0], numS*8, MPI_DOUBLE, IDS, tag, mpi_comm_level1, &request);
    }

    /* Recv */
    if (Num_Rcv_Grid1[IDR]!=0){

      numR = Num_Rcv_Grid1[IDR]*Ngrid2*Ngrid3;        

      tmp_array1 = (double*)malloc(sizeof(double)*numR*8);

      MPI_Recv(&tmp_array1[0], numR*8, MPI_DOUBLE, IDR, tag, mpi_comm_level1, &stat);

      for (i=0; i<Num_Rcv_Grid1[IDR]; i++){ 
      n1 = Rcv_Grid1[IDR][i];
        nn1 = My_Cell0[n1];
        nn0 = n1 - Start_Grid1[myid];
        MN0 = i*Ngrid2*Ngrid3;
        MN1 = nn0*Ngrid2*Ngrid3;
        for (n2=0; n2<Ngrid2; n2++){
          MN2 = n2*Ngrid3;
          for (n3=0; n3<Ngrid3; n3++){
            MNN1 = MN1 + MN2 + n3;
            MNN0 = MN0 + MN2 + n3;
            Den[0][MNN1] = tmp_array1[         MNN0];
            Den[1][MNN1] = tmp_array1[  numR + MNN0];
            ADen[MNN1]   = tmp_array1[2*numR + MNN0];
            PDen[MNN1]   = tmp_array1[3*numR + MNN0];
            dVHart[MNN1] = tmp_array1[4*numR + MNN0];
            Vxc[0][MNN1] = tmp_array1[5*numR + MNN0];
            Vxc[1][MNN1] = tmp_array1[6*numR + MNN0];
            Vef[MNN1]    = tmp_array1[7*numR + MNN0];
        }
      }
      }

      free(tmp_array1);
    }

    if (Num_Snd_Grid1[IDS]!=0){
      MPI_Wait(&request,&stat);
      free(tmp_array0);
    }
  }

  /* use own densities */

  for (n1=0; n1<Ngrid1; n1++){

    if (Cell_ID0[n1]==myid){

      nn0 = My_Cell0[n1];
      nn1 = n1 - Start_Grid1[myid];

      if (nn0!=-1 && Start_Grid1[myid]<=n1 && n1<=End_Grid1[myid]){

        MN0 = nn0*Ngrid2*Ngrid3;  
        MN1 = nn1*Ngrid2*Ngrid3;

        for (n2=0; n2<Ngrid2; n2++){
          MN2 = n2*Ngrid3;
          for (n3=0; n3<Ngrid3; n3++){
            MNN0 = MN0 + MN2 + n3;
            MNN1 = MN1 + MN2 + n3;

            Den[0][MNN1] = Density_Grid[0][MNN0];
            Den[1][MNN1] = Density_Grid[1][MNN0];
            ADen[MNN1]   = ADensity_Grid[MNN0];
            PDen[MNN1]   = PCCDensity_Grid[MNN0];
            dVHart[MNN1] = dVHart_Grid[MNN0];
            Vxc[0][MNN1] = Vxc_Grid[0][MNN0];
            Vxc[1][MNN1] = Vxc_Grid[1][MNN0];
            Vef[MNN1]    = VEF_Grid[MNN0];
          }    
        }
      }
    }
  }

  /****************************************************
           electric energy by electric field 
  ****************************************************/

  if (E_Field_switch==1){

    sum = 0.0;
    for (MN=0; MN<My_GNum; MN++){
      sum += (Den[0][MN] + Den[1][MN])*Vef[MN];
    }
    My_Eef = sum*GridVol;

    MPI_Barrier(mpi_comm_level1);
    MPI_Allreduce(&My_Eef, &Eef, 1, MPI_DOUBLE, MPI_SUM, mpi_comm_level1);

    ECE[12] = Eef;
  }
  else {
    ECE[12] = 0.0;
  }

  if (F_VEF_flag==0){
    ECE[12] = 0.0;
  }

  /****************************************************
     EH1 = 1/2\int \delta n(r) \delta V_H dr
  ****************************************************/

  sum = 0.0;
  for (MN=0; MN<My_GNum; MN++){
    sum += (Den[0][MN] + Den[1][MN] - ADen[MN])*dVHart[MN];
  }
  My_EH1 = 0.5*sum*GridVol;

  /*
  My_charge = 0.0;
  for (MN=0; MN<My_GNum; MN++){
    My_charge += Den[0][MN] + Den[1][MN];
  }
  My_charge = My_charge*GridVol;
  MPI_Allreduce(&My_charge, &sum_charge, 1, MPI_DOUBLE, MPI_SUM, mpi_comm_level1);
  printf("sum_charge=%15.12f\n",sum_charge);
  */

  /****************************************************
     EXC = \sum_{\sigma} n_{\sigma}\epsilon_{xc}
  ****************************************************/

  if      (SpinP_switch==0) max_spin = 0;
  else if (SpinP_switch==1) max_spin = 1;
  else if (SpinP_switch==3) max_spin = 1;

  if (PCC_switch==0){
    for (spin=0; spin<=max_spin; spin++){
      sum = 0.0;
      for (MN=0; MN<My_GNum; MN++){
        if (1.0e-14<Den[spin][MN]){
          sum += Den[spin][MN]*Vxc[spin][MN];
      }
      }
      My_EXC[spin] = sum*GridVol;
    }
  }
  else if (PCC_switch==1){
    for (spin=0; spin<=max_spin; spin++){
      sum = 0.0;
      for (MN=0; MN<My_GNum; MN++){
        if (1.0e-14<Den[spin][MN]){
          sum += (Den[spin][MN] + PDen[MN])*Vxc[spin][MN];
      }
      }
      My_EXC[spin] = sum*GridVol;
    }
  }

  /****************************************************
   MPI:

   EH1, EXC
  ****************************************************/

  MPI_Barrier(mpi_comm_level1);
  MPI_Allreduce(&My_EH1, &EH1, 1, MPI_DOUBLE, MPI_SUM, mpi_comm_level1);
  for (spin=0; spin<=max_spin; spin++){
    MPI_Allreduce(&My_EXC[spin], &EXC[spin], 1, MPI_DOUBLE, MPI_SUM, mpi_comm_level1);
  }

  if (SpinP_switch==0){
    ECE[5] = EH1;
    ECE[6] = EXC[0];
    ECE[7] = EXC[0];
  }
  else if (SpinP_switch==1 || SpinP_switch==3) {
    ECE[5] = EH1;
    ECE[6] = EXC[0];
    ECE[7] = EXC[1];
  }

  if (F_dVHart_flag==0){
    ECE[5] = 0.0;
  }

  if (F_Vxc_flag==0){
    ECE[6] = 0.0;
    ECE[7] = 0.0;
  }

  /****************************************************
           calculation of dipole moment
  ****************************************************/

  /* contribution from electron density */

  My_E_dpx = 0.0;
  My_E_dpy = 0.0;
  My_E_dpz = 0.0;

  My_E_dpx_BG = 0.0;
  My_E_dpy_BG = 0.0;
  My_E_dpz_BG = 0.0;

  for (MN=0; MN<My_GNum; MN++){

    n1 = MN/(Ngrid2*Ngrid3);
    n2 = (MN - n1*(Ngrid2*Ngrid3))/Ngrid3;
    n3 = MN - n1*(Ngrid2*Ngrid3) - n2*Ngrid3;
    n1 = n1 + Start_Grid1[myid];

    x = (double)n1*gtv[1][1] + (double)n2*gtv[2][1]
      + (double)n3*gtv[3][1] + Grid_Origin[1];
    y = (double)n1*gtv[1][2] + (double)n2*gtv[2][2]
      + (double)n3*gtv[3][2] + Grid_Origin[2];
    z = (double)n1*gtv[1][3] + (double)n2*gtv[2][3]
      + (double)n3*gtv[3][3] + Grid_Origin[3];

    den = Den[0][MN] + Den[1][MN];
       
    My_E_dpx += den*x;
    My_E_dpy += den*y;
    My_E_dpz += den*z; 

    My_E_dpx_BG += x;
    My_E_dpy_BG += y;
    My_E_dpz_BG += z; 
  }

  MPI_Allreduce(&My_E_dpx, &E_dpx, 1, MPI_DOUBLE, MPI_SUM, mpi_comm_level1);
  MPI_Allreduce(&My_E_dpy, &E_dpy, 1, MPI_DOUBLE, MPI_SUM, mpi_comm_level1);
  MPI_Allreduce(&My_E_dpz, &E_dpz, 1, MPI_DOUBLE, MPI_SUM, mpi_comm_level1);

  MPI_Allreduce(&My_E_dpx_BG, &E_dpx_BG, 1, MPI_DOUBLE, MPI_SUM, mpi_comm_level1);
  MPI_Allreduce(&My_E_dpy_BG, &E_dpy_BG, 1, MPI_DOUBLE, MPI_SUM, mpi_comm_level1);
  MPI_Allreduce(&My_E_dpz_BG, &E_dpz_BG, 1, MPI_DOUBLE, MPI_SUM, mpi_comm_level1);

  E_dpx = E_dpx*GridVol;
  E_dpy = E_dpy*GridVol;
  E_dpz = E_dpz*GridVol;

  cden_BG = system_charge/Cell_Volume; 

  E_dpx_BG = E_dpx_BG*GridVol*cden_BG;
  E_dpy_BG = E_dpy_BG*GridVol*cden_BG;
  E_dpz_BG = E_dpz_BG*GridVol*cden_BG;

  /* contribution from core charge */

  My_C_dpx = 0.0;
  My_C_dpy = 0.0;
  My_C_dpz = 0.0;

  for (Mc_AN=1; Mc_AN<=Matomnum; Mc_AN++){
    Gc_AN = M2G[Mc_AN];
    x = Gxyz[Gc_AN][1];
    y = Gxyz[Gc_AN][2];
    z = Gxyz[Gc_AN][3];

    spe = WhatSpecies[Gc_AN];
    charge = Spe_Core_Charge[spe];
    My_C_dpx += charge*x;
    My_C_dpy += charge*y;
    My_C_dpz += charge*z;
  }

  MPI_Allreduce(&My_C_dpx, &C_dpx, 1, MPI_DOUBLE, MPI_SUM, mpi_comm_level1);
  MPI_Allreduce(&My_C_dpy, &C_dpy, 1, MPI_DOUBLE, MPI_SUM, mpi_comm_level1);
  MPI_Allreduce(&My_C_dpz, &C_dpz, 1, MPI_DOUBLE, MPI_SUM, mpi_comm_level1);

  AU2Debye = 2.54174776;

  dipole_moment[0][1] = AU2Debye*(C_dpx - E_dpx - E_dpx_BG);
  dipole_moment[0][2] = AU2Debye*(C_dpy - E_dpy - E_dpy_BG);
  dipole_moment[0][3] = AU2Debye*(C_dpz - E_dpz - E_dpz_BG);

  dipole_moment[1][1] = AU2Debye*C_dpx;
  dipole_moment[1][2] = AU2Debye*C_dpy;
  dipole_moment[1][3] = AU2Debye*C_dpz;

  dipole_moment[2][1] = -AU2Debye*E_dpx;
  dipole_moment[2][2] = -AU2Debye*E_dpy;
  dipole_moment[2][3] = -AU2Debye*E_dpz;

  dipole_moment[3][1] = -AU2Debye*E_dpx_BG;
  dipole_moment[3][2] = -AU2Debye*E_dpy_BG;
  dipole_moment[3][3] = -AU2Debye*E_dpz_BG;

  AbsD = sqrt( dipole_moment[0][1]*dipole_moment[0][1]
             + dipole_moment[0][2]*dipole_moment[0][2]
             + dipole_moment[0][3]*dipole_moment[0][3] );

  if (myid==Host_ID){

    printf("\n*******************************************************\n"); fflush(stdout);
    printf("                  Dipole moment (Debye)                 \n");  fflush(stdout);
    printf("*******************************************************\n\n"); fflush(stdout);

    printf(" Absolute D %17.8f\n\n",AbsD);
    printf("                      Dx                Dy                Dz\n"); fflush(stdout);
    printf(" Total       %17.8f %17.8f %17.8f\n",
              dipole_moment[0][1],dipole_moment[0][2],dipole_moment[0][3]);fflush(stdout);
    printf(" Core        %17.8f %17.8f %17.8f\n",
              dipole_moment[1][1],dipole_moment[1][2],dipole_moment[1][3]);fflush(stdout);
    printf(" Electron    %17.8f %17.8f %17.8f\n",
              dipole_moment[2][1],dipole_moment[2][2],dipole_moment[2][3]);fflush(stdout);
    printf(" Back ground %17.8f %17.8f %17.8f\n",
              dipole_moment[3][1],dipole_moment[3][2],dipole_moment[3][3]);fflush(stdout);

    /********************************************************
             write the dipole moments to a file
    ********************************************************/

    fnjoint(filepath,filename,file_DPM);

    if ((fp_DPM = fopen(file_DPM,"w")) != NULL){

#ifdef xt3
      setvbuf(fp_DPM,buf,_IOFBF,fp_bsize);  /* setvbuf */
#endif

      fprintf(fp_DPM,"\n");
      fprintf(fp_DPM,"***********************************************************\n");
      fprintf(fp_DPM,"***********************************************************\n");
      fprintf(fp_DPM,"                    Dipole moment (Debye)                  \n");
      fprintf(fp_DPM,"***********************************************************\n");
      fprintf(fp_DPM,"***********************************************************\n\n");

      fprintf(fp_DPM," Absolute D %17.8f\n\n",AbsD);
      fprintf(fp_DPM,"                      Dx                Dy                Dz\n");
      fprintf(fp_DPM," Total       %17.8f %17.8f %17.8f\n",
                dipole_moment[0][1],dipole_moment[0][2],dipole_moment[0][3]);
      fprintf(fp_DPM," Core        %17.8f %17.8f %17.8f\n",
                dipole_moment[1][1],dipole_moment[1][2],dipole_moment[1][3]);
      fprintf(fp_DPM," Electron    %17.8f %17.8f %17.8f\n",
                dipole_moment[2][1],dipole_moment[2][2],dipole_moment[2][3]);
      fprintf(fp_DPM," Back ground %17.8f %17.8f %17.8f\n",
                dipole_moment[3][1],dipole_moment[3][2],dipole_moment[3][3]);

      fclose(fp_DPM);
    }
    else{
      printf("Failure of saving the DPM file.\n");fflush(stdout);
    }
  }

  /****************************************************
   freeing of arrays:

    free(Den[0]);
    free(Den[1]);
    free(Den);
    free(ADen);
    free(PDen);
    free(dVHart);
    free(Vxc[0]);
    free(Vxc[1]);
    free(Vxc);
    free(Vef);
  ****************************************************/

  free(Den[0]);
  free(Den[1]);
  free(Den);
  free(ADen);
  free(PDen);
  free(dVHart);
  free(Vxc[0]);
  free(Vxc[1]);
  free(Vxc);
  free(Vef);
}








void EH0_TwoCenter(int Gc_AN, int h_AN, double VH0ij[4])
{ 
  int n1,ban;
  int Gh_AN,Rn,wan1,wan2;
  double dv,x,y,z,r,r2,va0,rho0,dr_va0;
  double z2,sum,sumr,sumx,sumy,sumz,wt;

  Gh_AN = natn[Gc_AN][h_AN];
  Rn = ncn[Gc_AN][h_AN];
  wan1 = WhatSpecies[Gc_AN];
  ban = Spe_Spe2Ban[wan1];
  wan2 = WhatSpecies[Gh_AN];
  dv = dv_EH0[ban];
  
  sum = 0.0;
  sumr = 0.0;

  for (n1=0; n1<TGN_EH0[ban]; n1++){
    x = GridX_EH0[ban][n1];
    y = GridY_EH0[ban][n1];
    z = GridZ_EH0[ban][n1];
    rho0 = Arho_EH0[ban][n1];
    wt = Wt_EH0[ban][n1];
    z2 = z - Dis[Gc_AN][h_AN];
    r2 = sqrt(x*x + y*y + z2*z2);

    /* for empty atoms or finite elemens basis */
    if (r2<1.0e-10) r2 = 1.0e-10;

    va0 = VH_AtomF(wan2,r2);
    sum = sum + wt*va0*rho0;

    if (h_AN!=0 && 1.0e-14<r2){
      dr_va0 = Dr_VH_AtomF(wan2,r2);
      sumr = sumr - wt*rho0*dr_va0*z2/r2;
    }
  }

  sum  = sum*dv;

  if (h_AN!=0){

    /* for empty atoms or finite elemens basis */
    r = Dis[Gc_AN][h_AN];
    if (r<1.0e-10) r = 1.0e-10;

    x = Gxyz[Gc_AN][1] - (Gxyz[Gh_AN][1] + atv[Rn][1]);
    y = Gxyz[Gc_AN][2] - (Gxyz[Gh_AN][2] + atv[Rn][2]);
    z = Gxyz[Gc_AN][3] - (Gxyz[Gh_AN][3] + atv[Rn][3]);
    sumr = sumr*dv;
    sumx = sumr*x/r;
    sumy = sumr*y/r;
    sumz = sumr*z/r;
  }
  else{
    sumx = 0.0;
    sumy = 0.0;
    sumz = 0.0;
  }

  VH0ij[0] = sum;
  VH0ij[1] = sumx;
  VH0ij[2] = sumy;
  VH0ij[3] = sumz;
}





void EH0_TwoCenter_at_Cutoff(int wan1, int wan2, double VH0ij[4])
{ 
  int n1,ban;
  double dv,x,y,z,r1,r2,va0,rho0,dr_va0,rcut;
  double z2,sum,sumr,sumx,sumy,sumz,wt;

  ban = Spe_Spe2Ban[wan1];
  dv = dv_EH0[ban];
  rcut = Spe_Atom_Cut1[wan1] + Spe_Atom_Cut1[wan2];  

  sum = 0.0;
  for (n1=0; n1<TGN_EH0[ban]; n1++){
    x = GridX_EH0[ban][n1];
    y = GridY_EH0[ban][n1];
    z = GridZ_EH0[ban][n1];
    rho0 = Arho_EH0[ban][n1];
    wt = Wt_EH0[ban][n1];
    z2 = z - rcut;
    r2 = sqrt(x*x + y*y + z2*z2);
    va0 = VH_AtomF(wan2,r2);
    sum = sum + wt*va0*rho0;
  }

  sum  = sum*dv;
  sumx = 0.0;
  sumy = 0.0;
  sumz = 0.0;

  VH0ij[0] = sum;
  VH0ij[1] = sumx;
  VH0ij[2] = sumy;
  VH0ij[3] = sumz;
}










double Calc_Ehub()
{   
 /****************************************************
         LDA+U energy correction added by MJ
  ****************************************************/

  int Mc_AN,Gc_AN,wan1;
  int cnt1,cnt2,l1,mul1,m1,l2,mul2,m2;
  int spin,max_spin;
  double My_Ehub,Ehub,Uvalue;
  int numprocs,myid,ID;

  /* MPI */
  MPI_Comm_size(mpi_comm_level1,&numprocs);
  MPI_Comm_rank(mpi_comm_level1,&myid);

 /****************************************************
                 caculation of My_Ehub
  ****************************************************/

  if      (SpinP_switch==0) max_spin = 0;
  else if (SpinP_switch==1) max_spin = 1;
  else if (SpinP_switch==3) max_spin = 1;

  My_Ehub = 0.0;

  for (Mc_AN=1; Mc_AN<=Matomnum; Mc_AN++){
    Gc_AN = M2G[Mc_AN];
    wan1 = WhatSpecies[Gc_AN];

   /****************************************************
                     collinear case
    ****************************************************/

    if (SpinP_switch!=3){

      for (spin=0; spin<=max_spin; spin++){


        /* Hubbard term, 0.5*Tr(N) */

      cnt1 = 0;
      for(l1=0; l1<=Spe_MaxL_Basis[wan1]; l1++ ){
        for(mul1=0; mul1<Spe_Num_Basis[wan1][l1]; mul1++){

          Uvalue = Hub_U_Basis[wan1][l1][mul1];
          for(m1=0; m1<(2*l1+1); m1++){

            My_Ehub += 0.5*Uvalue*DM_onsite[0][spin][Mc_AN][cnt1][cnt1];

            cnt1++;
          }
        }
      }


        /* Hubbard term, -0.5*Tr(N*N) */

      cnt1 = 0;
      for(l1=0; l1<=Spe_MaxL_Basis[wan1]; l1++ ){
        for(mul1=0; mul1<Spe_Num_Basis[wan1][l1]; mul1++){
          for(m1=0; m1<(2*l1+1); m1++){

            cnt2 = 0;
            for(l2=0; l2<=Spe_MaxL_Basis[wan1]; l2++ ){
            for(mul2=0; mul2<Spe_Num_Basis[wan1][l2]; mul2++){
              for(m2=0; m2<(2*l2+1); m2++){

                if (l1==l2 && mul1==mul2){
                  Uvalue = Hub_U_Basis[wan1][l1][mul1];
                  My_Ehub -= 0.5*Uvalue*DM_onsite[0][spin][Mc_AN][cnt1][cnt2]*
                                      DM_onsite[0][spin][Mc_AN][cnt2][cnt1];
                }

                cnt2++;
              }
            }
            }

            cnt1++;
          }
        }
      }
      }
    }

   /****************************************************
                     non-collinear case
    ****************************************************/

    else {

      /* Hubbard term, 0.5*Tr(N) */

      cnt1 = 0;
      for(l1=0; l1<=Spe_MaxL_Basis[wan1]; l1++ ){
      for(mul1=0; mul1<Spe_Num_Basis[wan1][l1]; mul1++){

        Uvalue = Hub_U_Basis[wan1][l1][mul1];
        for(m1=0; m1<(2*l1+1); m1++){

          My_Ehub += 0.5*Uvalue*( NC_OcpN[0][0][0][Mc_AN][cnt1][cnt1].r
                          + NC_OcpN[0][1][1][Mc_AN][cnt1][cnt1].r);

          cnt1++;
        }
      }
      }

      /* Hubbard term, -0.5*Tr(N*N) */

      cnt1 = 0;
      for(l1=0; l1<=Spe_MaxL_Basis[wan1]; l1++ ){
      for(mul1=0; mul1<Spe_Num_Basis[wan1][l1]; mul1++){
        for(m1=0; m1<(2*l1+1); m1++){

          cnt2 = 0;
          for(l2=0; l2<=Spe_MaxL_Basis[wan1]; l2++ ){
            for(mul2=0; mul2<Spe_Num_Basis[wan1][l2]; mul2++){
            for(m2=0; m2<(2*l2+1); m2++){

              if (l1==l2 && mul1==mul2){

                Uvalue = Hub_U_Basis[wan1][l1][mul1];

                My_Ehub -= 0.5*Uvalue*( NC_OcpN[0][0][0][Mc_AN][cnt1][cnt2].r*
                                  NC_OcpN[0][0][0][Mc_AN][cnt1][cnt2].r
                                  +
                                  NC_OcpN[0][0][0][Mc_AN][cnt1][cnt2].i*
                                  NC_OcpN[0][0][0][Mc_AN][cnt1][cnt2].i
                                  +
                                  NC_OcpN[0][0][1][Mc_AN][cnt1][cnt2].r*
                                  NC_OcpN[0][0][1][Mc_AN][cnt1][cnt2].r
                                  +
                                  NC_OcpN[0][0][1][Mc_AN][cnt1][cnt2].i*
                                  NC_OcpN[0][0][1][Mc_AN][cnt1][cnt2].i
                                  +
                                  NC_OcpN[0][1][0][Mc_AN][cnt1][cnt2].r*
                                  NC_OcpN[0][1][0][Mc_AN][cnt1][cnt2].r
                                  +
                                  NC_OcpN[0][1][0][Mc_AN][cnt1][cnt2].i*
                                  NC_OcpN[0][1][0][Mc_AN][cnt1][cnt2].i
                                  +
                                  NC_OcpN[0][1][1][Mc_AN][cnt1][cnt2].r*
                                  NC_OcpN[0][1][1][Mc_AN][cnt1][cnt2].r
                                  +
                                  NC_OcpN[0][1][1][Mc_AN][cnt1][cnt2].i*
                                  NC_OcpN[0][1][1][Mc_AN][cnt1][cnt2].i );

              }

              cnt2++;
            }
            }
          }

          cnt1++;
        }
      }
      }
    }
  }

  if (SpinP_switch==0) My_Ehub = 2.0*My_Ehub;

 /****************************************************
                      MPI My_Ehub
  ****************************************************/

  MPI_Allreduce(&My_Ehub, &Ehub, 1, MPI_DOUBLE, MPI_SUM, mpi_comm_level1);

  /* if (F_U_flag==0) */
  if (F_U_flag==0) Ehub = 0.0;
  
  return Ehub;  
}

Generated by  Doxygen 1.6.0   Back to index