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

Make_ESM.c

/**********************************************************************
  Make_ESM.c:

     Make_ESM.c is a subroutine to make an exchange splitting matrix
     which is used for evaluation of exchange coupling constants.
     This calculation is valid for collinear spin polarized case.

  Log of Make_ESM.c:

     22/Jan/2004  Released by T.Ozaki

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

#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

void Make_ESM()
{
  double time0;
  int Mc_AN,Gc_AN,Mh_AN,h_AN,Gh_AN;
  int q_AN,Gq_AN,Mq_AN,Rnq,Qwan,Nq;
  int i,j,k,Cwan,Hwan,NO0,NO1;
  int Rnh,Rnk,spin,N,NumC[4];
  int n1,n2,n3,L0,Mul0,M0,L1,Mul1,M1;
  int Nc,GNc,GRc,Nog,Nh,MN,XC_P_switch;
  int calc_on,ESM_Lc_AN,tno0,tno1,num;
  int *ESM_size;
  double x,y,z,dx,dy,dz,fw;
  double bc,dv,r,theta,phi,sum,tmp0;
  double xo,yo,zo,S_coordinate[3];
  double Cxyz[4];
  double **ChiV0;
  double *FuzzyW;
  double TStime,TEtime;
  double *tmp_array0;
  double *tmp_array1;
  int numprocs,myid,tag=999,ID;

  MPI_Status stat;
  MPI_Request request;

  MPI_Comm_size(mpi_comm_level1,&numprocs);
  MPI_Comm_rank(mpi_comm_level1,&myid);

  dtime(&TStime);

  if (myid==Host_ID) printf("<Make_ESM>  calculate exchange splitting matrices\n");

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

    double ChiV0[List_YOUSO[7]][List_YOUSO[11]];
    double FuzzyW[List_YOUSO[11]];
    int ESM_size[Num_ESM];
  ****************************************************/

  ChiV0 = (double**)malloc(sizeof(double*)*List_YOUSO[7]);
  for (j=0; j<List_YOUSO[7]; j++){
    ChiV0[j] = (double*)malloc(sizeof(double)*List_YOUSO[11]);
  }

  FuzzyW = (double*)malloc(sizeof(double)*List_YOUSO[11]);
  ESM_size = (int*)malloc(sizeof(int)*Num_ESM);

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

   double ESM[Num_ESM][FNAN][FNAN][tno0][tno1]
  ****************************************************/

  ESM = (double*****)malloc(sizeof(double****)*Num_ESM);
  for (k=0; k<Num_ESM; k++){

    ESM_size[k] = 0;

    Gc_AN = ESM_Atoms[k];
    ESM[k] = (double****)malloc(sizeof(double***)*(FNAN[Gc_AN]+1));

    for (h_AN=0; h_AN<=FNAN[Gc_AN]; h_AN++){
      Gh_AN = natn[Gc_AN][h_AN];
      Hwan = WhatSpecies[Gh_AN];
      tno0 = Spe_Total_CNO[Hwan];
      ESM[k][h_AN] = (double***)malloc(sizeof(double**)*(FNAN[Gc_AN]+1));

      for (q_AN=0; q_AN<=FNAN[Gc_AN]; q_AN++){
        Gq_AN = natn[Gc_AN][q_AN];
        Qwan = WhatSpecies[Gq_AN];
        tno1 = Spe_Total_CNO[Qwan];
        ESM[k][h_AN][q_AN] = (double**)malloc(sizeof(double*)*tno0);

        for (i=0; i<tno0; i++){
          ESM[k][h_AN][q_AN][i] = (double*)malloc(sizeof(double)*tno1);
      }
        
        ESM_size[k] += tno0*tno1;
      }
    }
  }

  /*
  if (Cnt_kind==1) {
    Contract_Hamiltonian( H, CntH, OLP, CntOLP );
    if (SO_switch==1) Contract_iHNL(iHNL,iCntHNL);
  }
  */

  /*****************************************************
    calculation of matrix elements for dVH + Vxc + VNA
  *****************************************************/

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

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

    /* check a list of ESM_Atoms */    
    calc_on = 0;  
    k = 0;
    do {
      if (Gc_AN==ESM_Atoms[k]){
        calc_on = 1;
        ESM_Lc_AN = k;
      }
      k++; 
    } while(calc_on==0 && k<Num_ESM);

    if (calc_on==1){ 

      /* calculate fuzzy weight */
      for (Nog=0; Nog<NumOLG[Mc_AN][0]; Nog++){
      GNc = GListTAtoms0[Mc_AN][0][Nog];
      GRc = GListTCells0[Mc_AN][0][Nog];
      Get_Grid_XYZ(GNc,Cxyz);
      x = Cxyz[1] + atv[GRc][1];
      y = Cxyz[2] + atv[GRc][2]; 
      z = Cxyz[3] + atv[GRc][3];
      Nc = GListTAtoms1[Mc_AN][0][Nog];
        FuzzyW[Nc] = Fuzzy_Weight(Gc_AN,0,x,y,z);
      }

      /* for the first orbital */
      for (h_AN=0; h_AN<=FNAN[Gc_AN]; h_AN++){

      Gh_AN = natn[Gc_AN][h_AN];
      Mh_AN = F_G2M[Gh_AN];
      Rnh = ncn[Gc_AN][h_AN];
      Hwan = WhatSpecies[Gh_AN];
      NO0 = Spe_Total_CNO[Hwan];

      /* multiply orbitals of Mh_AN by splitting potential on site Mc_AN */
      for (i=0; i<NO0; i++){
        for (Nc=0; Nc<GridN_Atom[Gc_AN]; Nc++){
          ChiV0[i][Nc] = 0.0;
        }
      }

        for (Nog=0; Nog<NumOLG[Mc_AN][h_AN]; Nog++){
        Nc = GListTAtoms1[Mc_AN][h_AN][Nog];
        Nh = GListTAtoms2[Mc_AN][h_AN][Nog];
        MN = GListTAtoms3[Mc_AN][h_AN][Nog];

        for (i=0; i<NO0; i++){


          ChiV0[i][Nc] = Orbs_Grid[Mh_AN][i][Nh]*FuzzyW[Nc]*
                           (Vpot_Grid[0][MN]-Vpot_Grid[1][MN]);

          /*
          ChiV0[i][Nc] = Orbs_Grid[Mh_AN][i][Nh]*
                           (Vpot_Grid[0][MN]-Vpot_Grid[1][MN]);
          */

        }
      }

      /* for the second orbital */
      for (q_AN=0; q_AN<=FNAN[Gc_AN]; q_AN++){

        Gq_AN = natn[Gc_AN][q_AN];
        Mq_AN = F_G2M[Gq_AN];
        Rnq = ncn[Gc_AN][q_AN];
        Qwan = WhatSpecies[Gq_AN];
        NO1 = Spe_Total_CNO[Qwan];

        /*
        printf("Mc_AN h_AN q_AN %2d %2d %2d\n",Mc_AN,h_AN,q_AN);
        */

        for (i=0; i<NO0; i++){
          for (j=0; j<NO1; j++){

            sum = 0.0;
            for (Nog=0; Nog<NumOLG[Mc_AN][q_AN]; Nog++){
            Nc = GListTAtoms1[Mc_AN][q_AN][Nog];
            Nq = GListTAtoms2[Mc_AN][q_AN][Nog];
            sum += ChiV0[i][Nc]*Orbs_Grid[Mq_AN][j][Nq];
            }

            /*
            printf("%10.6f ",sum*GridVol);
            */

              ESM[ESM_Lc_AN][h_AN][q_AN][i][j] = sum*GridVol;

          }

          /*
          printf("\n");
          */

        }

      }
      }
    }
  }

  /****************************************************
    MPI ESM
  ****************************************************/

  for (k=0; k<Num_ESM; k++){
    Gc_AN = ESM_Atoms[k];

    tmp_array0 = (double*)malloc(sizeof(double)*ESM_size[k]);
    tmp_array1 = (double*)malloc(sizeof(double)*ESM_size[k]);

    if (G2ID[Gc_AN]!=Host_ID){

      /* send */

      num = -1; 
      for (h_AN=0; h_AN<=FNAN[Gc_AN]; h_AN++){

      Gh_AN = natn[Gc_AN][h_AN];
      Mh_AN = F_G2M[Gh_AN];
      Rnh = ncn[Gc_AN][h_AN];
      Hwan = WhatSpecies[Gh_AN];
      NO0 = Spe_Total_CNO[Hwan];

      for (q_AN=0; q_AN<=FNAN[Gc_AN]; q_AN++){

        Gq_AN = natn[Gc_AN][q_AN];
        Mq_AN = F_G2M[Gq_AN];
        Rnq = ncn[Gc_AN][q_AN];
        Qwan = WhatSpecies[Gq_AN];
        NO1 = Spe_Total_CNO[Qwan];

        for (i=0; i<NO0; i++){
          for (j=0; j<NO1; j++){
              num++;
              tmp_array0[num] = ESM[k][h_AN][q_AN][i][j];
          }
        }
      }
      }      

      /* MPI_Isend */
      MPI_Isend(&tmp_array0[0], ESM_size[k], MPI_DOUBLE, Host_ID,
                tag, mpi_comm_level1, &request);
      MPI_Wait(&request,&stat);

      /* receive */

      if (myid==Host_ID){

      /* MPI_Recv */
      MPI_Recv(&tmp_array1[0], ESM_size[k], MPI_DOUBLE,
             G2ID[Gc_AN], tag, mpi_comm_level1, &stat);

      num = -1; 
      for (h_AN=0; h_AN<=FNAN[Gc_AN]; h_AN++){

        Gh_AN = natn[Gc_AN][h_AN];
        Mh_AN = F_G2M[Gh_AN];
        Rnh = ncn[Gc_AN][h_AN];
        Hwan = WhatSpecies[Gh_AN];
        NO0 = Spe_Total_CNO[Hwan];

        for (q_AN=0; q_AN<=FNAN[Gc_AN]; q_AN++){

          Gq_AN = natn[Gc_AN][q_AN];
          Mq_AN = F_G2M[Gq_AN];
          Rnq = ncn[Gc_AN][q_AN];
          Qwan = WhatSpecies[Gq_AN];
          NO1 = Spe_Total_CNO[Qwan];

          for (i=0; i<NO0; i++){
            for (j=0; j<NO1; j++){
            num++;
            ESM[k][h_AN][q_AN][i][j] = tmp_array1[num];
            }
          }
        }
      }
      }
      
    }
   
    free(tmp_array0);
    free(tmp_array1);
  }

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

    double ChiV0[List_YOUSO[7]][List_YOUSO[11]];
    int ESM_size[Num_ESM];
  ****************************************************/

  for (j=0; j<List_YOUSO[7]; j++){
    free(ChiV0[j]);
  }
  free(ChiV0);

  free(FuzzyW);
  free(ESM_size);

  /* for time */
  dtime(&TEtime);
  time0 = TEtime - TStime;

}


Generated by  Doxygen 1.6.0   Back to index