abacus-develop/source/source_estate/module_charge/symmetry_rho.cpp at c365a3afd18cced4dfdc424c120fe0a0e3bfb4a6 · deepmodeling/abacus-develop · GitHub

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
#include "symmetry_rho.h"

#include "source_hamilt/module_xc/xc_functional.h"

Symmetry_rho::Symmetry_rho()
{
}

Symmetry_rho::~Symmetry_rho()
{
}

void Symmetry_rho::symmetrize_rho(const int nspin,
                                   const Charge& chr,
                                   const ModulePW::PW_Basis* pw,
                                   ModuleSymmetry::Symmetry& symm)
{
    Symmetry_rho srho;
    for (int is = 0; is < nspin; is++)
    {
        srho.begin(is, chr, pw, symm);
    }
}

void Symmetry_rho::begin(const int& spin_now,
                         const Charge& chr,
                         const ModulePW::PW_Basis* rho_basis,
                         ModuleSymmetry::Symmetry& symm) const
{
    assert(spin_now < 4); // added by zhengdy-soc

	if (ModuleSymmetry::Symmetry::symm_flag != 1)
	{
		return;
	}

    ModuleBase::TITLE("Symmetry_rho", "begin");
    ModuleBase::timer::tick("Symmetry_rho","begin");

// both parallel and serial
// if(symm.nrot==symm.nrotk) //pure point-group, do rho_symm in real space
// {
// 	psymm(chr.rho[spin_now], rho_basis, Pgrid, symm);
// 	if(XC_Functional::get_ked_flag()) psymm(chr.kin_r[spin_now],
// rho_basis,Pgrid,symm);
// }
// else	//space group, do rho_symm in reciprocal space

    rho_basis->real2recip(chr.rho[spin_now], chr.rhog[spin_now]);

    psymmg(chr.rhog[spin_now], rho_basis, symm); // need to modify

    rho_basis->recip2real(chr.rhog[spin_now], chr.rho[spin_now]);

    if (XC_Functional::get_ked_flag() || chr.cal_elf)
    {
        // Use std::vector to manage kin_g instead of raw pointer
        std::vector<std::complex<double>> kin_g(chr.ngmc);
        rho_basis->real2recip(chr.kin_r[spin_now], kin_g.data());
        psymmg(kin_g.data(), rho_basis, symm);
        rho_basis->recip2real(kin_g.data(), chr.kin_r[spin_now]);
    }

    ModuleBase::timer::tick("Symmetry_rho","begin");
    return;
}

void Symmetry_rho::begin(const int& spin_now,
                         double** rho,
                         std::complex<double>** rhog,
                         int ngmc,
                         double** kin_r,
                         const ModulePW::PW_Basis* rho_basis,
                         ModuleSymmetry::Symmetry& symm) const
{
    assert(spin_now < 4); // added by zhengdy-soc

    if (ModuleSymmetry::Symmetry::symm_flag != 1)
    {
        return;
    }

    ModuleBase::TITLE("Symmetry_rho", "begin");
    ModuleBase::timer::tick("Symmetry_rho","begin");

    // both parallel and serial
    // if(symm.nrot==symm.nrotk) //pure point-group, do rho_symm in real space
    // {
    // 	psymm(CHR.rho[spin_now], rho_basis, Pgrid, symm);
    // 	if(XC_Functional::get_ked_flag()) psymm(CHR.kin_r[spin_now],
    // rho_basis,Pgrid,symm);
    // }
    // else	//space group, do rho_symm in reciprocal space
    {
        rho_basis->real2recip(rho[spin_now], rhog[spin_now]);
        psymmg(rhog[spin_now], rho_basis, symm);
        rho_basis->recip2real(rhog[spin_now], rho[spin_now]);

        if (XC_Functional::get_ked_flag())
        {
            // Use std::vector to manage kin_g instead of raw pointer
            std::vector<std::complex<double>> kin_g(ngmc);
            rho_basis->real2recip(kin_r[spin_now], kin_g.data());
            psymmg(kin_g.data(), rho_basis, symm);
            rho_basis->recip2real(kin_g.data(), kin_r[spin_now]);
        }
    }

    ModuleBase::timer::tick("Symmetry_rho","begin");
    return;
}

void Symmetry_rho::psymm(double* rho_part,
                         const ModulePW::PW_Basis* rho_basis,
                         Parallel_Grid& Pgrid,
                         ModuleSymmetry::Symmetry& symm) const
{
    ModuleBase::TITLE("Symmetry_rho", "psymm");
    ModuleBase::timer::tick("Symmetry_rho","psymm");

#ifdef __MPI
    // reduce all rho from the first pool.
    std::vector<double> rhotot;
    if (GlobalV::MY_RANK == 0)
    {
        rhotot.resize(rho_basis->nxyz);
        ModuleBase::GlobalFunc::ZEROS(rhotot.data(), rho_basis->nxyz);
    }
    Pgrid.reduce(rhotot.data(), rho_part, false);

    if (GlobalV::MY_RANK == 0)
    {
        symm.rho_symmetry(rhotot.data(), rho_basis->nx, rho_basis->ny, rho_basis->nz);
#else
    symm.rho_symmetry(rho_part, rho_basis->nx, rho_basis->ny, rho_basis->nz);
#endif
        /*
        int count = 0;
        GlobalV::ofs_running << scientific;
        for(int iz=0; iz<rho_basis->nz; iz++)
        {
            GlobalV::ofs_running << "\n iz=" << iz;
            for(int iy=0; iy<rho_basis->ny; iy++)
            {
                for(int ix=0; ix<rho_basis->nx; ix++)
                {
                    if(count%5==0) GlobalV::ofs_running << "\n";
                    ++count;
                    GlobalV::ofs_running << " " << rhotot[ix*rho_basis->ny*rho_basis->nz+iy*rho_basis->nz+iz];
                }
            }
        }
        */
#ifdef __MPI
    }

	Pgrid.bcast(rhotot.data(), rho_part, GlobalV::MY_RANK);
#endif

    ModuleBase::timer::tick("Symmetry_rho","psymm");
    return;
}