Skip to content

Commit

Permalink
Implemented the tadpole diagram
Browse files Browse the repository at this point in the history
  • Loading branch information
@ttadano
ttadano committed Feb 23, 2016
1 parent e0cac06 commit d2c7f1c
Show file tree
Hide file tree
Showing 3 changed files with 79 additions and 4 deletions.
17 changes: 13 additions & 4 deletions anphon/relaxation.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -113,6 +113,8 @@ void Relaxation::setup_relaxation()
if (calc_realpart && integration->ismear != 0) {
error->exit("setup_relaxation", "Sorry. REALPART = 1 can be used only with ISMEAR = 0");
}

dynamical->modify_eigenvectors();
}

if (kpoint->kpoint_mode == 2) {
Expand Down Expand Up @@ -1232,6 +1234,7 @@ void Relaxation::perform_mode_analysis()

double *damping_a;
double omega_shift;
std::complex<double> *self_tadpole;
std::complex<double> *self_a, *self_b, *self_c, *self_d, *self_e;
std::complex<double> *self_f, *self_g, *self_h, *self_i, *self_j;

Expand Down Expand Up @@ -1262,6 +1265,7 @@ void Relaxation::perform_mode_analysis()
memory->allocate(damping_a, NT);
memory->allocate(self_a, NT);
memory->allocate(self_b, NT);
memory->allocate(self_tadpole, NT);

if (quartic_mode == 2) {
memory->allocate(self_c, NT);
Expand Down Expand Up @@ -1353,6 +1357,7 @@ void Relaxation::perform_mode_analysis()

if (calc_realpart) {

selfenergy->selfenergy_tadpole(NT, T_arr, omega, knum, snum, self_tadpole);
selfenergy->selfenergy_a(NT, T_arr, omega, knum, snum, self_a);

if (quartic_mode == 1) {
Expand All @@ -1373,22 +1378,26 @@ void Relaxation::perform_mode_analysis()
ofs_shift << std::endl;
ofs_shift << "# mode = " << snum + 1<< std::endl;
ofs_shift << "# Frequency = " << writes->in_kayser(omega) << std::endl;
ofs_shift << "## T[K], Shift3 (cm^-1) (bubble)";
ofs_shift << "## T[K], Shift3 (cm^-1) (tadpole), Shift3 (cm^-1) (bubble)";
if (quartic_mode == 1) ofs_shift << ", Shift4 (cm^-1) (loop)";
ofs_shift << ", Shifted frequency (cm^-1)";
ofs_shift << std::endl;


for (j = 0; j < NT; ++j) {
ofs_shift << std::setw(10) << T_arr[j] << std::setw(15) << writes->in_kayser(-self_a[j].real());
ofs_shift << std::setw(10) << T_arr[j];
ofs_shift << std::setw(15) << writes->in_kayser(-self_tadpole[j].real());
ofs_shift << std::setw(15) << writes->in_kayser(-self_a[j].real());

omega_shift = omega - self_a[j].real();
omega_shift = omega - self_tadpole[j].real() - self_a[j].real();

if (quartic_mode == 1) {
ofs_shift << std::setw(15) << writes->in_kayser(-self_b[j].real());
omega_shift -= self_b[j].real();
}
ofs_shift << std::setw(15) << writes->in_kayser(omega_shift);
ofs_shift << std::endl;

}

ofs_shift.close();
Expand All @@ -1400,6 +1409,7 @@ void Relaxation::perform_mode_analysis()
memory->deallocate(damping_a);
memory->deallocate(self_a);
memory->deallocate(self_b);
memory->deallocate(self_tadpole);

if (quartic_mode == 2) {
memory->deallocate(self_c);
Expand Down Expand Up @@ -2482,7 +2492,6 @@ void Relaxation::setup_quartic()
}
}

dynamical->modify_eigenvectors();

memory->deallocate(invsqrt_mass_p);
}
Expand Down
65 changes: 65 additions & 0 deletions anphon/selfenergy.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -65,6 +65,71 @@ void Selfenergy::mpi_reduce_complex(unsigned int N, std::complex<double> *in_mpi
#endif
}

void Selfenergy::selfenergy_tadpole(const unsigned int N, double *T, const double omega, const unsigned int knum, const unsigned int snum, std::complex<double> *ret)
{
unsigned int i;
unsigned int ik2;
unsigned int is1, is2;
unsigned int arr_cubic1[3], arr_cubic2[3];
std::complex<double> v3_tmp1, v3_tmp2;
std::complex<double> *ret_mpi, *ret_tmp;

double T_tmp;
double n2;
double omega1, omega2;
double factor;

arr_cubic1[0] = ns * kpoint->knum_minus[knum] + snum;
arr_cubic1[1] = ns * knum + snum;

memory->allocate(ret_mpi, N);
memory->allocate(ret_tmp, N);

for (i = 0; i < N; ++i) ret[i] = std::complex<double>(0.0, 0.0);

for (is1 = 0; is1 < ns; ++is1) {
arr_cubic1[2] = is1;
arr_cubic2[0] = is1;
omega1 = dynamical->eval_phonon[0][is1];

if (omega1 < eps8) continue;

v3_tmp1 = relaxation->V3(arr_cubic1);

for (i = 0; i < N; ++i) ret_mpi[i] = std::complex<double>(0.0, 0.0);

for (ik2 = mympi->my_rank; ik2 < nk; ik2 += mympi->nprocs) {
for (is2 = 0; is2 < ns; ++is2) {
arr_cubic2[1] = ns * ik2 + is2;
arr_cubic2[2] = ns * kpoint->knum_minus[ik2] + is2;

v3_tmp2 = relaxation->V3(arr_cubic2);
omega2 = dynamical->eval_phonon[ik2][is2];

if (omega2 < eps8) continue;

for (i = 0; i < N; ++i) {
T_tmp = T[i];
n2 = thermodynamics->fB(omega2, T_tmp);
// ret[i] += v3_tmp1 * v3_tmp2 * (2.0*n2 + 1.0) / omega1;
ret_mpi[i] += v3_tmp2 * (2.0 * n2 + 1.0);
}
}
}
mpi_reduce_complex(N, ret_mpi, ret_tmp);

for (i = 0; i < N; ++i) {
ret[i] += ret_tmp[i] * v3_tmp1 / omega1;
}
}

factor = -1.0 / (static_cast<double>(nk) * std::pow(2.0, 3));
for (i = 0; i < N; ++i) ret[i] *= factor;

memory->deallocate(ret_tmp);
memory->deallocate(ret_mpi);
}

void Selfenergy::selfenergy_a(const unsigned int N, double *T, const double omega, const unsigned int knum, const unsigned int snum, std::complex<double> *ret)
{
/*
Expand Down
1 change: 1 addition & 0 deletions anphon/selfenergy.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -22,6 +22,7 @@ namespace PHON_NS {
~Selfenergy();

void setup_selfenergy();
void selfenergy_tadpole(const unsigned int, double *, const double, const unsigned int, const unsigned int, std::complex<double> *);
void selfenergy_a(const unsigned int, double *, const double, const unsigned int, const unsigned int, std::complex<double> *);
void selfenergy_b(const unsigned int, double *, const double, const unsigned int, const unsigned int, std::complex<double> *);
void selfenergy_c(const unsigned int, double *, const double, const unsigned int, const unsigned int, std::complex<double> *);
Expand Down

0 comments on commit d2c7f1c

Please sign in to comment.