Generated from mg_eig_state.c with ROBODoc v3.2.2 on Fri Sep 14 14:23:19 2001
TABLE OF CONTENTS
- QMD-MGDFT/mg_eig_state.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 mg_eig_state(STATE *sp, int tid, REAL *vtot)
there are two routines, one for real (Gamma point only) and the other for complex
Multigrid solver for a single electronic orbital.
Algorithm:
1. Eigenvalues are computed once per cycle using the Mehrstallen
operator.
2. Presmoothing of wavefunctions on the finest grid level using
the Mehrstallen operator.
3. The defect from the last Mehrstallen smoothing is then used
as the right hand side for a richardson defect correction
scheme using a 7-point discretization of the del-squared operator.
This is multigridded.
4. Mehrstallen post smoothing.
INPUTS
sp: point to orbital structure STATE (see md.h)
tid: thread ID
vtot: total potential = vxc +vnuc +vh
OUTPUT
wave functions and eigenvalues are updated and stored in STATE
PARENTS
scf.c threads.c
CHILDREN
gather_psi.c app_nl.c pack_ptos.c app_cilr.c genvpsi.c app_cir.c global_sums.c
trade_images.c app_smooth.c mgrid_solv.c pack_stop_axpy.c scatter_psi.c
SOURCE
#include <float.h>
#include <math.h>
#include <stdlib.h>
#include <stdio.h>
#include "md.h"
#if GAMMA_PT
void mg_eig_state(STATE *sp, int tid, REAL *vtot)
{
int idx, cycles;
int nits, pbasis, sbasis;
REAL eig, diag, t1, t2, t3, t4;
REAL *work1, *work2;
REAL *tmp_psi, *nvtot, *res, *rhs, *sg_psi, *sg_twovpsi;
int eig_pre[6] = {0, 3, 6, 2, 2, 2};
int eig_post[6] = {0, 3, 6, 2, 2, 2};
int ione=1;
int dimx, dimy, dimz, levels;
REAL hxgrid, hygrid, hzgrid, sb_step;
REAL tarr[8];
#if 1
REAL time1, time2;
time1 = my_crtc();
#endif
nits = ct.eig_parm.gl_pre + ct.eig_parm.gl_pst;
dimx = sp->dimx;
dimy = sp->dimy;
dimz = sp->dimz;
hxgrid = sp->hxgrid;
hygrid = sp->hygrid;
hzgrid = sp->hzgrid;
levels = ct.eig_parm.levels;
sb_step = 1.0;
pbasis = sp->pbasis;
sbasis = sp->sbasis;
/* Grab some memory */
sg_psi = (REAL *)get_mem(12 * (sbasis));
tmp_psi = sg_psi + sbasis;
res = tmp_psi + sbasis;
work2 = res + 2 * sbasis;
sg_twovpsi = work2 + 4*sbasis;
nvtot = sg_twovpsi + sbasis;
rhs = nvtot + sbasis;
work1 = rhs + sbasis;
gather_psi(tmp_psi, NULL, sp, tid);
/* Get the non-local operator acting on psi and store in nvtot */
app_nl(tmp_psi, NULL, nvtot, NULL, sp->istate, sp->firstflag, sp->kidx, tid);
/* Smoothing cycles */
for(cycles = 0;cycles <= nits;cycles++) {
/* Pack psi into smoothing array */
pack_ptos(sg_psi, tmp_psi, dimx, dimy, dimz);
/* Apply Mehrstellen left and right hand operators */
diag = app_cilr(sg_psi, work2, res, dimx, dimy, dimz,
hxgrid, hygrid, hzgrid);
diag = -1.0 / diag;
/* Generate 2 * V * psi */
genvpsi(tmp_psi, sg_twovpsi, vtot, nvtot, NULL, 0.0,
dimx, dimy, dimz);
/* B operating on 2*V*psi stored in work1 */
app_cir(sg_twovpsi, work1, dimx, dimy, dimz);
t1 = -ONE;
QMD_saxpy(pbasis, t1, work2, ione, work1, ione);
/* If this is the first time through compute the eigenvalue */
if(cycles == 0) {
eig = 0.0;
t2 = 0.0;
for(idx = 0;idx < pbasis;idx++){
t2 += tmp_psi[idx] * res[idx];
eig += tmp_psi[idx] * work1[idx];
} /* end for */
tarr[0] = t2;
tarr[1] = eig;
idx = 2;
global_sums(tarr, &idx);
eig = tarr[1] / (TWO * tarr[0]);
sp->eig = eig;
} /* end if */
/* Now either smooth the wavefunction or do a multigrid cycle */
if(cycles == ct.eig_parm.gl_pre) {
t1 = TWO * eig;
for(idx = 0;idx < P0_BASIS;idx++) {
res[idx] = t1 * res[idx] - work1[idx];
} /* end for */
/* Pack the residual data into multigrid array */
pack_ptos(sg_psi, res, dimx, dimy, dimz);
trade_images(sg_psi, dimx, dimy, dimz, pct.neighbors);
/* Smooth it once and store the smoothed residual in work1 */
t1 = 145.0;
app_smooth((S0_GRID *)sg_psi, (S0_GRID *)work1, t1);
/* Do multigrid step with solution returned in sg_twovpsi */
mgrid_solv(sg_twovpsi, work1, work2,
dimx, dimy, dimz, hxgrid,
hygrid, hzgrid,
0, pct.neighbors, levels,
eig_pre, eig_post, 1, sb_step);
/* The correction is in a smoothing grid so we use this
* routine to update the orbital which is stored in a physical grid.
*/
t1 = -1.0;
pack_stop_axpy(sg_twovpsi, tmp_psi, t1, dimx, dimy, dimz);
}
else {
t1 = TWO * eig;
t2 = ZERO;
t4 = ct.eig_parm.gl_step * diag;
if(cycles == 0)t4 = 2.2 * diag;
if(cycles == 1)t4 = 1.1 * diag;
for(idx = 0;idx < P0_BASIS;idx++) {
t3 = t1 * res[idx] - work1[idx];
t2 += t3 * t3;
tmp_psi[idx] += t4 * t3;
} /* end for */
if(cycles == 0) {
t2 = real_sum_all(t2);
t1 = (REAL)(ct.psi_nbasis);
sp->res = ct.hmaxgrid * ct.hmaxgrid * sqrt(t2 / t1) * 0.25;
} /* end if */
} /* end if */
} /* end for */
/* Copy psi back into regular array */
scatter_psi(tmp_psi, NULL, sp, tid);
sp->firstflag++;
/* Release our memory */
release_mem(sg_psi);
#if 1
time2 = my_crtc();
md_timings(MG_EIGTIME, (time2 - time1), tid);
#endif
} /* end mg_eig_state */
#else
/* Complex version */
void mg_eig_state(STATE *sp, int tid, REAL *vtot)
{
int idx, cycles;
int nits, pbasis, sbasis;
REAL eig=0.0, eigR, eigI, diag, t2, t2R, t2I;
REAL t1, *kdr, *kdi, *work1R, *work2R, *work1I, *work2I, *gx, *gy, *gz;
REAL *tmp_psiR, *nvtotR, *resR, *sg_psiR, *sg_twovpsiR;
REAL *tmp_psiI, *nvtotI, *resI, *sg_psiI, *sg_twovpsiI;
int eig_pre[6] = {0, 2, 2, 2, 2, 2};
int eig_post[6] = {0, 2, 2, 2, 2, 2};
int ione=1;
int dimx, dimy, dimz, levels;
REAL hxgrid, hygrid, hzgrid, sb_step;
#if 1
REAL time1, time2;
time1 = my_crtc();
#endif
nits = ct.eig_parm.gl_pre + ct.eig_parm.gl_pst;
dimx = sp->dimx;
dimy = sp->dimy;
dimz = sp->dimz;
hxgrid = sp->hxgrid;
hygrid = sp->hygrid;
hzgrid = sp->hzgrid;
levels = ct.eig_parm.levels;
sb_step = 1.0;
pbasis = sp->pbasis;
sbasis = sp->sbasis;
/* Grab some memory */
sg_psiR = (REAL *)get_mem(12 * (sbasis));
tmp_psiR = sg_psiR + sbasis;
resR = tmp_psiR + sbasis;
work2R = resR + 2 * sbasis;
sg_twovpsiR = work2R + 4*sbasis;
nvtotR = sg_twovpsiR + sbasis;
work1R = nvtotR + sbasis;
sg_psiI = (REAL *)get_mem(17 * (sbasis));
tmp_psiI = sg_psiI + sbasis;
resI = tmp_psiI + sbasis;
work2I = resI + 2 * sbasis;
sg_twovpsiI = work2I + 4*sbasis;
nvtotI = sg_twovpsiI + sbasis;
work1I = nvtotI + sbasis;
gx = work1I + sbasis;
gy = gx + sbasis;
gz = gy + sbasis;
kdr = gz + sbasis;
kdi = kdr + sbasis;
gather_psi(tmp_psiR, tmp_psiI, sp, tid);
/* Get the non-local operator acting on psi and store in nvtot */
app_nl(tmp_psiR, tmp_psiI, nvtotR, nvtotI, sp->istate, sp->firstflag, sp->kidx, tid);
/* Smoothing cycles */
for(cycles = 0;cycles <= nits;cycles++) {
/* First do the real part of the orbitals */
/* Pack psi into smoothing arrays */
pack_ptos(sg_psiR, tmp_psiR, dimx, dimy, dimz);
/* Apply Mehrstellen left and right hand operators */
diag = app_cilr(sg_psiR, work2R, resR, dimx, dimy, dimz,
hxgrid, hygrid, hzgrid);
pack_ptos(sg_psiI, tmp_psiI, dimx, dimy, dimz);
/* Apply Mehrstellen left and right hand operators */
diag = app_cilr(sg_psiI, work2I, resI, dimx, dimy, dimz,
hxgrid, hygrid, hzgrid);
diag = -1.0 / diag;
/* Apply the gradient operator to psi */
app_grad((S0_GRID *)sg_psiI, (P0_GRID *)gx, (P0_GRID *)gy, (P0_GRID *)gz);
for(idx = 0;idx < P0_BASIS;idx++) {
kdr[idx] = (ct.kp[sp->kidx].kvec[0] * gx[idx] +
ct.kp[sp->kidx].kvec[1] * gy[idx] +
ct.kp[sp->kidx].kvec[2] * gz[idx]);
} /* end for */
/* Generate 2 * V * psi */
genvpsi(tmp_psiR, sg_twovpsiR, vtot, nvtotR, kdr, ct.kp[sp->kidx].kmag,
dimx, dimy, dimz);
/* B operating on 2*V*psi stored in work1 */
app_cir(sg_twovpsiR, work1R, dimx, dimy, dimz);
t1 = -ONE;
QMD_saxpy(pbasis, t1, work2R, ione, work1R, ione);
/* Apply the gradient operator to psi */
app_grad((S0_GRID *)sg_psiR, (P0_GRID *)gx, (P0_GRID *)gy, (P0_GRID *)gz);
for(idx = 0;idx < P0_BASIS;idx++) {
kdi[idx] = -(ct.kp[sp->kidx].kvec[0] * gx[idx] +
ct.kp[sp->kidx].kvec[1] * gy[idx] +
ct.kp[sp->kidx].kvec[2] * gz[idx]);
} /* end for */
/* Generate 2 * V * psi */
genvpsi(tmp_psiI, sg_twovpsiI, vtot, nvtotI, kdi, ct.kp[sp->kidx].kmag,
dimx, dimy, dimz);
/* B operating on 2*V*psi stored in work1 */
app_cir(sg_twovpsiI, work1I, dimx, dimy, dimz);
t1 = -ONE;
QMD_saxpy(pbasis, t1, work2I, ione, work1I, ione);
/* If this is the first time through compute the eigenvalue */
if(cycles == 0) {
eigR = t2R = ZERO;
for(idx = 0;idx < pbasis;idx++){
t2R += tmp_psiR[idx] * resR[idx];
eigR += tmp_psiR[idx] * work1R[idx];
} /* end for */
eigI = t2I = ZERO;
for(idx = 0;idx < pbasis;idx++){
t2I += tmp_psiI[idx] * resI[idx];
eigI += tmp_psiI[idx] * work1I[idx];
} /* end for */
eigR = real_sum_all(eigR);
t2R = real_sum_all(t2R);
eigI = real_sum_all(eigI);
t2I = real_sum_all(t2I);
sp->eig = (eigR + eigI) / (TWO * (t2R + t2I));
eig = sp->eig;
} /* end if */
/* Next we have to generate the residual vector for smoothing */
/* or multigridding. */
t1 = TWO * eig;
for(idx = 0;idx < P0_BASIS;idx++) {
resR[idx] = t1 * resR[idx] - work1R[idx];
resI[idx] = t1 * resI[idx] - work1I[idx];
} /* end for */
if(cycles == 0) {
t1 = snrm2(&pbasis, resR, &ione);
t1 = t1 * t1;
t1 = real_sum_all(t1);
t2 = (REAL)(ct.psi_nbasis);
sp->res = ct.hmaxgrid * ct.hmaxgrid * sqrt(t1 / t2) * 0.25;
} /* end if */
/* Now either smooth the wavefunction or do a multigrid cycle */
if(cycles == ct.eig_parm.gl_pre) {
/* Real components */
/* Pack the residual data into multigrid array */
pack_ptos(sg_psiR, resR, dimx, dimy, dimz);
trade_images(sg_psiR, dimx, dimy, dimz, pct.neighbors);
/* Smooth it once and store the smoothed residual in work1 */
t1 = 145.0;
app_smooth((S0_GRID *)sg_psiR, (S0_GRID *)work1R, t1);
/* Do multigrid step with solution returned in sg_twovpsi */
mgrid_solv(sg_twovpsiR, work1R, work2R,
dimx, dimy, dimz, hxgrid,
hygrid, hzgrid,
0, pct.neighbors, levels,
eig_pre, eig_post, 1, sb_step);
/* The correction is in a smoothing grid so we use this
* routine to update the orbital which is stored in a physical grid.
*/
t1 = -1.0;
pack_stop_axpy(sg_twovpsiR, tmp_psiR, t1, dimx, dimy, dimz);
/* Imaginary components */
/* Pack the residual data into multigrid array */
pack_ptos(sg_psiI, resI, dimx, dimy, dimz);
trade_images(sg_psiI, dimx, dimy, dimz, pct.neighbors);
/* Smooth it once and store the smoothed residual in work1 */
t1 = 145.0;
app_smooth((S0_GRID *)sg_psiI, (S0_GRID *)work1I, t1);
/* Do multigrid step with solution returned in sg_twovpsi */
mgrid_solv(sg_twovpsiI, work1I, work2I,
dimx, dimy, dimz, hxgrid,
hygrid, hzgrid,
0, pct.neighbors, levels,
eig_pre, eig_post, 1, sb_step);
/* The correction is in a smoothing grid so we use this
* routine to update the orbital which is stored in a physical grid.
*/
t1 = -1.1;
pack_stop_axpy(sg_twovpsiI, tmp_psiI, t1, dimx, dimy, dimz);
}
else {
t1 = ct.eig_parm.gl_step * diag;
if(cycles == 0)t1 = 2.2 * diag;
if(cycles == 1)t1 = 1.1 * diag;
/* Update wavefuntion */
QMD_saxpy(pbasis, t1, resR, ione, tmp_psiR, ione);
QMD_saxpy(pbasis, t1, resI, ione, tmp_psiI, ione);
} /* end if */
} /* end for */
/* Copy psi back into regular array */
scatter_psi(tmp_psiR, tmp_psiI, sp, tid);
sp->firstflag++;
/* Release our memory */
release_mem(sg_psiI);
release_mem(sg_psiR);
#if 1
time2 = my_crtc();
md_timings(MG_EIGTIME, (time2 - time1), tid);
#endif
} /* end mg_eig_state */
#endif