Generated from get_rho.c with ROBODoc v3.2.2 on Fri Sep 14 14:23:17 2001
TABLE OF CONTENTS
- QMD-MGDFT/get_rho.c
NAME
Ab initio real space code with multigrid acceleration
Quantum molecular dynamics package.
Version: 2.1.5
COPYRIGHT
Copyright (C) 1995 Emil Briggs
Copyright (C) 1998 Emil Briggs, Charles Brabec, Mark Wensell,
Dan Sullivan, Chris Rapcewicz, Jerzy Bernholc
Copyright (C) 2001 Emil Briggs, Wenchang Lu,
Marco Buongiorno Nardelli,Charles Brabec,
Mark Wensell,Dan Sullivan, Chris Rapcewicz,
Jerzy Bernholc
FUNCTION
void get_rho(STATE *states, REAL *rho)
Generates new charge density and mix linearly with old one.
INPUTS
states: point to orbital structure (see md.h)
rho: old charge density
OUTPUT
rho: updated charge density
PARENTS
scf.c
CHILDREN
gather_psi.c symmetrize_rho.f
SOURCE
/*
get_rho.c
*/
#include <float.h>
#include <math.h>
#include <stdlib.h>
#include <stdio.h>
#include "md.h"
#if (MACHINE_TYPE != PARALLEL_SMP)
void get_rho(STATE *states, REAL *rho)
{
int istate, kpt, n, incx, idx;
STATE *sp;
REAL t1, *work, *tmp_psi;
work = get_mem(P0_BASIS);
tmp_psi = get_mem(P0_BASIS);
if(work == NULL)
error_handler("get_rho", "Out of memory");
/* scale charge accumulator */
n = P0_BASIS;
incx = 1;
t1 = ONE - ct.mix;
QMD_sscal(n, t1, rho, incx);
for(kpt = 0;kpt < ct.num_kpts;kpt++) {
sp = ct.kp[kpt].kstate;
/* Loop over states and accumulate charge */
for(istate = 0;istate < ct.num_states;istate++) {
n = P0_BASIS;
incx = 1;
t1 = ct.mix * sp->occupation * ct.kp[kpt].kweight;
gather_psi(tmp_psi, NULL, sp, 0);
for(idx = 0;idx < P0_BASIS;idx++) {
work[idx] = tmp_psi[idx] * tmp_psi[idx];
} /* end for */
QMD_saxpy(n, t1, work, incx, rho, incx);
#if !GAMMA_PT
gather_psi(NULL, tmp_psi, sp, 0);
for(idx = 0;idx < P0_BASIS;idx++) {
work[idx] = tmp_psi[idx] * tmp_psi[idx];
} /* end for */
QMD_saxpy(n, t1, work, incx, rho, incx);
#endif
sp++;
} /* end for */
} /* end for */
#if !GAMMA_PT
/* Symmetrize the density */
symmetrize_rho((P0_GRID *)rho);
#endif
/* Check total charge. */
ct.tcharge = ZERO;
for(idx = 0;idx < P0_BASIS;idx++) {
ct.tcharge += rho[idx];
} /* end for */
ct.tcharge = real_sum_all(ct.tcharge);
ct.tcharge = ct.tcharge * ct.vel;
/* Renormalize charge */
t1 = ct.nel / ct.tcharge;
for(idx = 0;idx < P0_BASIS;idx++) {
rho[idx] = rho[idx] * t1;
}
/* release our memory */
release_mem(tmp_psi);
release_mem(work);
} /* end get_rho */
#else
/* Must be a parallel SMP or NUMA machine. We assume that this is
* always called directly after an orthogonalization step (the ortho
* routines store thread local copies of the density which are used here.
*/
void get_rho(STATE *states, REAL *rho)
{
int idx, thread, offset, ione=1, stop;
REAL *rptr;
ct.tcharge = ZERO;
for(thread = 0;thread < ct.num_threads;thread++) {
ct.tcharge += thread_control[thread].tcharge;
}
ct.tcharge = ct.vel * ct.tcharge;
for(idx = 0;idx < P0_BASIS;idx++) rho[idx] = 0.0;
offset = 0;
for(thread = 0;thread < ct.num_threads;thread++) {
stop = thread_control[thread].numpt;
rptr = thread_control[thread].rho;
QMD_scopy( stop, rptr, ione, &rho[offset + ct.node_space_offset], ione);
offset += stop;
} /* end for */
#if !GAMMA_PT
/* Symmetrize the density */
symmetrize_rho((P0_GRID *)rho);
#endif
} /* end get_rho */
void get_rho_smp(SCF_THREAD_CONTROL *s)
{
int kpt, koffset, istate, n, incx, idx;
STATE *sp;
REAL t1, *work, *tmp_psiR, *tmp_psiI;
work = get_mem(2 * s->numpt);
if(work == NULL)
error_handler("get_rho", "Out of memory");
/* scale charge accumulator */
n = s->numpt;
incx = 1;
t1 = ONE - ct.mix;
QMD_sscal(n, t1, s->rho, incx);
for(kpt = 0;kpt < ct.num_kpts;kpt++) {
/* Loop over states and accumulate charge */
koffset = 2 * kpt * s->lda * ct.num_states;
tmp_psiR = &thread_control[s->tid].base_mem[koffset];
for(istate = 0;istate < ct.num_states;istate++) {
sp = &s->my_states[kpt*ct.num_states + istate];
t1 = ct.mix * sp->occupation * ct.kp[kpt].kweight;
#if GAMMA_PT
for(idx = 0;idx < s->numpt;idx++) {
work[idx] = tmp_psiR[idx] * tmp_psiR[idx];
} /* end for */
QMD_saxpy(n, t1, work, incx, s->rho, incx);
tmp_psiR += s->lda;
#else
tmp_psiI = tmp_psiR + s->lda;
for(idx = 0;idx < s->numpt;idx++) {
work[idx] = tmp_psiR[idx] * tmp_psiR[idx] +
tmp_psiI[idx] * tmp_psiI[idx];
} /* end for */
QMD_saxpy(n, t1, work, incx, s->rho, incx);
tmp_psiR += 2*s->lda;
#endif
} /* end for */
} /* end for */
/* Check total charge. */
s->tcharge = ZERO;
for(idx = 0;idx < s->numpt;idx++) {
s->tcharge += s->rho[idx];
} /* end for */
/* release our memory */
release_mem(work);
} /* end get_rho_smp */
#endif