Generated from lforce.c with ROBODoc v3.2.2 on Fri Sep 14 14:23:19 2001
TABLE OF CONTENTS
- QMD-MGDFT/lforce.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 lforce(REAL *rho, REAL *vh)
Evaluates the ionic force component due to the local part of potentials
and the compensating charges via the Hellman-Feynman theorem.
INPUTS
rho: total charge density
vh: Hartree potential
OUTPUT
forces are added to structure ct.ions
PARENTS
force.c
CHILDREN
get_index.c to_cartisian.c linint.c
SOURCE
#include <float.h>
#include <math.h>
#include <stdio.h>
#include <stdlib.h>
#include "md.h"
#include "inline.h"
void lforce(REAL *rho, REAL *vh)
{
int ix, iy, iz;
int ion, idx;
int *pvec, docount;
int ilow, jlow, klow, ihi, jhi, khi, map;
int icut, itmp, icenter;
int Aix[NX_GRID], Aiy[NY_GRID], Aiz[NZ_GRID];
REAL r, xc, yc, zc, Zv, rc, rcnorm, t1;
REAL ax[3], bx[3], x, y, z, rc2, invdr, norm1;
REAL fx, fy, fz;
SPECIES *sp;
ION *iptr;
REAL *rx, *ry, *rz;
REAL *urx, *ury, *urz;
#if 1
REAL time1, time2;
time1 = my_crtc();
#endif
rx = get_mem(7 * P0_BASIS);
ry = rx + P0_BASIS;
rz = ry + P0_BASIS;
urx = rz + P0_BASIS;
ury = urx + P0_BASIS;
urz = ury + P0_BASIS;
pvec = (int *)get_mem(P0_BASIS);
/* 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];
icenter = sp->ldim / 2;
icut = (icenter + 1) * (icenter + 1);
Zv = sp->zvalence;
rc = sp->rc;
rc2 = rc * rc;
rcnorm = rc*rc*rc * pow(PI, 1.5);
rcnorm = ONE / rcnorm;
rc2 = sp->rc * sp->rc;
norm1 = - Zv * rcnorm / rc2;
invdr = ONE / sp->drlig;
/* 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 */
my_barrier();
docount = 0;
if(map) {
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)) ) {
/* Cut it off if required */
itmp = (ix - icenter)*(ix - icenter) +
(iy - icenter)*(iy - icenter) +
(iz - icenter)*(iz - icenter);
if(icut >= itmp) {
pvec[docount] =
PY0_GRID * PZ0_GRID * (Aix[ix] % PX0_GRID) +
PZ0_GRID * (Aiy[iy] % PY0_GRID) +
(Aiz[iz] % PZ0_GRID);
ax[0] = xc - iptr->xtal[0];
ax[1] = yc - iptr->xtal[1];
ax[2] = zc - iptr->xtal[2];
to_cartesian(ax, bx);
x = bx[0];
y = bx[1];
z = bx[2];
r = metric(ax);
t1 = linint(sp->drlocalig, r, invdr);
r += 1.0e-10;
rx[docount] = -t1 * x / r;
ry[docount] = -t1 * y / r;
rz[docount] = -t1 * z / r;
t1 = 2.0 * norm1 * exp(-r*r / rc2);
urx[docount] = x * t1;
ury[docount] = y * t1;
urz[docount] = z * t1;
docount++;
} /* end if */
} /* end if */
zc += ct.hzgrid;
} /* end for */
yc += ct.hygrid;
} /* end for */
xc += ct.hxgrid;
} /* end for */
} /* end if */
my_barrier();
fx = fy = fz = 0.0;
for(idx = 0;idx < docount;idx++) {
fx += rx[idx] * rho[pvec[idx]];
fy += ry[idx] * rho[pvec[idx]];
fz += rz[idx] * rho[pvec[idx]];
fx += urx[idx] * vh[pvec[idx]];
fy += ury[idx] * vh[pvec[idx]];
fz += urz[idx] * vh[pvec[idx]];
} /* end for */
iptr->force[ct.fpt[0]][0] -= ct.vel * real_sum_all(fx);
iptr->force[ct.fpt[0]][1] -= ct.vel * real_sum_all(fy);
iptr->force[ct.fpt[0]][2] -= ct.vel * real_sum_all(fz);
} /* end for */
release_mem(pvec);
release_mem(rx);
#if 1
time2 = my_crtc();
md_timings(LFORCE_TIME, (time2 - time1), 0);
#endif
} /* end lforce */