Generated from init_nuc.c with ROBODoc v3.2.2 on Fri Sep 14 14:23:18 2001
TABLE OF CONTENTS
- QMD-MGDFT/init_nuc.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 init_nuc(REAL *vnuc, REAL *rhoc, REAL *rhocore)
Set up the nuclear local potentials, the core charges and the gaussian
compensating charges.
INPUTS
nothing
OUTPUT
vnuc: local part of pseudopotential
rhoc: gaussian compensating charges
rhocore: core charge density
PARENTS
cdfastrlx.c fastrlx.c init.c moldyn.c
CHILDREN
get_index.c
SOURCE
#include <float.h>
#include <math.h>
#include <stdio.h>
#include <stdlib.h>
#include "md.h"
#include "inline.h"
#include "efield.h"
void init_nuc(REAL *vnuc, REAL *rhoc, REAL *rhocore)
{
int ix, iy, iz;
int ion, idx;
int stop, ilow, jlow, klow, ihi, jhi, khi, map;
int Aix[NX_GRID], Aiy[NY_GRID], Aiz[NZ_GRID];
int thread, offset, icount;
int *pvec;
REAL r, xc, yc, zc, Zv, rc, rcnorm, t1;
REAL x[3], rc2, invdr;
SPECIES *sp;
ION *iptr;
/* Grab some memory for temporary storage */
pvec = get_mem( P0_BASIS );
/* Zero out the nuclear potential array */
for(idx = 0;idx < P0_BASIS;idx++) {
vnuc[idx] = ZERO;
} /* end for */
/* Initialize the compensating charge array and the core charge array */
for(idx = 0;idx < P0_BASIS;idx++) {
rhoc[idx] = ct.background_charge/P0_BASIS/ct.vel/NPES;
rhocore[idx] = ZERO;
} /* end for */
/* Loop over ions */
for(ion = 0;ion < ct.num_ions;ion++) {
/* Generate ion pointer */
iptr = &ct.ions[ion];
/* Get species type */
sp = ct.sp[iptr->species];
Zv = sp->zvalence;
rc = sp->rc;
rc2 = rc * rc;
rcnorm = rc*rc*rc * pow(PI, 1.5);
rcnorm = ONE / rcnorm;
/* Determine mapping indices or even if a mapping exists */
map = get_index(iptr, Aix, Aiy, Aiz, &ilow, &ihi, &jlow, &jhi, &klow, &khi,
sp->ldim, PX0_GRID, PY0_GRID, PZ0_GRID,
ct.psi_nxgrid, ct.psi_nygrid, ct.psi_nzgrid,
&iptr->lxcstart, &iptr->lycstart, &iptr->lzcstart);
/* If there is any overlap then we have to generate the mapping */
if(map) {
invdr = ONE / sp->drlig;
icount = 0;
xc = iptr->lxcstart;
for(ix = 0;ix < sp->ldim;ix++) {
yc = iptr->lycstart;
for(iy = 0;iy < sp->ldim;iy++) {
zc = iptr->lzcstart;
for(iz = 0;iz < sp->ldim;iz++) {
if( ((Aix[ix] >= ilow) && (Aix[ix] <= ihi)) &&
((Aiy[iy] >= jlow) && (Aiy[iy] <= jhi)) &&
((Aiz[iz] >= klow) && (Aiz[iz] <= khi)) ) {
pvec[icount] =
PY0_GRID * PZ0_GRID * (Aix[ix] % PX0_GRID) +
PZ0_GRID * (Aiy[iy] % PY0_GRID) +
(Aiz[iz] % PZ0_GRID);
x[0] = xc - iptr->xtal[0];
x[1] = yc - iptr->xtal[1];
x[2] = zc - iptr->xtal[2];
r = metric(x);
rhoc[pvec[icount]] += Zv * exp(-r*r / rc2) * rcnorm;
vnuc[pvec[icount]] +=
linint(&sp->localig[0], r, invdr);
if(sp->nlccflag)
rhocore[pvec[icount]] +=
linint(&sp->rhocorelig[0], r, invdr);
icount++;
} /* end if */
zc += ct.hzgrid;
} /* end for */
yc += ct.hygrid;
} /* end for */
xc += ct.hxgrid;
} /* end for */
} /* end if */
} /* end for */
/* Check compensating charges */
ct.crho = ZERO;
for(idx = 0;idx < P0_BASIS;idx++) {
ct.crho += rhoc[idx];
} /* end for */
ct.crho = ct.crho * ct.vel;
ct.crho = real_sum_all(ct.crho);
t1 = ZERO;
for(idx = 0;idx < P0_BASIS;idx++) {
if(rhocore[idx] < 0.0)rhocore[idx] = 0.0;
t1 += rhocore[idx];
} /* end for */
t1 = t1 * ct.vel;
t1 = real_sum_all(t1);
/* Release our memory */
release_mem(pvec);
/* Wait until everyone gets here */
my_barrier();
/* add the saw-tooth potential induced by applied electric field */
init_efield(vnuc);
#if (MACHINE_TYPE == PARALLEL_SMP)
offset = 0;
for(thread = 0;thread < ct.num_threads;thread++){
stop = thread_control[thread].numpt;
for(idx = 0;idx < stop;idx++) {
thread_control[thread].rhocore[idx] = rhocore[idx+offset];
}
offset += stop;
} /* end for */
#endif
} /* end init_nuc */