Generated from nlforce_s.c with ROBODoc v3.2.2 on Fri Sep 14 14:23:20 2001
TABLE OF CONTENTS
- QMD-MGDFT/nlforce_s.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 nlforce_d(int ion, int ip, int icount, STATE *states, int *dvec, int *pvec)
Calculates contribution to the ionic forces from a
set of s-projectors.
INPUTS
ion: ion index
ip: momentum l = 0 for s orbital
icount: number of grid point in the sphere
states: point to orbital structure (see md.h)
devc: true indicates the point is in the sphere, other not
pvec: index array of grid around an atom to physical grid
OUTPUT
forces are added in ct.ion structure
PARENTS
force.c
CHILDREN
metric.c to_cartesian.c linint.c gather_psi.c global_sum.c
SOURCE
#include <float.h>
#include <math.h>
#include <stdio.h>
#include <stdlib.h>
#include "md.h"
#include "inline.h"
/* Determines forces from an s-projector */
void nlforce_s(int ion, int ip, int icount, STATE *states, int *dvec, int *pvec)
{
int idx, ix, iy, iz, docount, st, incx = 1;
int stop, kpt;
REAL kweight;
REAL r, ax[3], bx[3], xc, yc, zc;
REAL t1, invdr;
REAL forcesR[4];
REAL forcesI[4];
REAL *tmp_psiR, *tmp_psiI, *loc_psiR, *loc_psiI, *prjptr;
REAL *rx, *ry, *rz;
REAL *pR, *pI;
SPECIES *sp;
ION *iptr;
STATE *stptr;
kpt = states->kidx;
kweight = ct.kp[kpt].kweight;
stop = pct.idxptrlen[ion];
pR = pct.phaseptr[ion];
pR += 2 * kpt * stop;
pI = pR + stop;
rx = get_mem(6 * P0_BASIS);
ry = rx + P0_BASIS;
rz = ry + P0_BASIS;
tmp_psiR = rz + P0_BASIS;
loc_psiR = tmp_psiR + P0_BASIS;
prjptr = loc_psiR + P0_BASIS;
#if !GAMMA_PT
tmp_psiI = get_mem(2*P0_BASIS);
loc_psiI = tmp_psiI + P0_BASIS;
#endif
/* Get ion pointer */
iptr = &ct.ions[ion];
/* Get species type */
sp = ct.sp[iptr->species];
invdr = ONE / sp->drnlig;
idx = docount = 0;
xc = iptr->nlxcstart;
for(ix = 0;ix < sp->nldim;ix++) {
yc = iptr->nlycstart;
for(iy = 0;iy < sp->nldim;iy++) {
zc = iptr->nlzcstart;
for(iz = 0;iz < sp->nldim;iz++) {
if(dvec[idx]) {
ax[0] = xc - iptr->xtal[0];
ax[1] = yc - iptr->xtal[1];
ax[2] = zc - iptr->xtal[2];
r = metric(ax);
to_cartesian(ax, bx);
t1 = linint(&sp->drkbplig[ip][0], r, invdr);
r += 1.0e-10;
rx[docount] = t1 * bx[0] / r;
ry[docount] = t1 * bx[1] / r;
rz[docount] = t1 * bx[2] / r;
prjptr[docount] = linint(&sp->kbplig[ip][0], r, invdr);
docount++;
} /* end if */
idx++;
zc += ct.hzgrid;
} /* end for */
yc += ct.hygrid;
} /* end for */
xc += ct.hxgrid;
} /* end for */
if(docount != icount)
error_handler("nlforce_s", "Problem with non-local force");
for(st = 0;st < ct.num_states;st++) {
stptr = &states[st];
/* Gather up the wavefunction */
#if GAMMA_PT
gather_psi(tmp_psiR, NULL, stptr, 0);
for(idx = 0;idx < icount;idx++) loc_psiR[idx] = tmp_psiR[pvec[idx]];
forcesR[0] = QMD_sdot(icount, rx, incx, loc_psiR, incx);
forcesR[1] = QMD_sdot(icount, ry, incx, loc_psiR, incx);
forcesR[2] = QMD_sdot(icount, rz, incx, loc_psiR, incx);
forcesR[3] = QMD_sdot(icount, prjptr, incx, loc_psiR, incx);
idx = 4;
global_sums(forcesR, &idx);
forcesR[0] = ct.vel * forcesR[0];
forcesR[1] = ct.vel * forcesR[1];
forcesR[2] = ct.vel * forcesR[2];
t1 = -2.0 * ct.vel * sp->kbnorm[ip] * stptr->occupation * forcesR[3];
iptr->force[ct.fpt[0]][0] -= t1 * forcesR[0];
iptr->force[ct.fpt[0]][1] -= t1 * forcesR[1];
iptr->force[ct.fpt[0]][2] -= t1 * forcesR[2];
#else
gather_psi(tmp_psiR, tmp_psiI, stptr, 0);
for(idx = 0;idx < icount;idx++) {
loc_psiR[idx] = tmp_psiR[pvec[idx]] * pR[idx] -
tmp_psiI[pvec[idx]] * pI[idx];
}
for(idx = 0;idx < icount;idx++) {
loc_psiI[idx] = tmp_psiI[pvec[idx]] * pR[idx] +
tmp_psiR[pvec[idx]] * pI[idx];
}
forcesR[0] = QMD_sdot(icount, rx, incx, loc_psiR, incx);
forcesR[1] = QMD_sdot(icount, ry, incx, loc_psiR, incx);
forcesR[2] = QMD_sdot(icount, rz, incx, loc_psiR, incx);
forcesR[3] = QMD_sdot(icount, prjptr, incx, loc_psiR, incx);
forcesI[0] = QMD_sdot(icount, rx, incx, loc_psiI, incx);
forcesI[1] = QMD_sdot(icount, ry, incx, loc_psiI, incx);
forcesI[2] = QMD_sdot(icount, rz, incx, loc_psiI, incx);
forcesI[3] = QMD_sdot(icount, prjptr, incx, loc_psiI, incx);
idx = 4;
global_sums(forcesR, &idx);
global_sums(forcesI, &idx);
forcesR[0] = ct.vel * forcesR[0];
forcesR[1] = ct.vel * forcesR[1];
forcesR[2] = ct.vel * forcesR[2];
forcesI[0] = ct.vel * forcesI[0];
forcesI[1] = ct.vel * forcesI[1];
forcesI[2] = ct.vel * forcesI[2];
t1 = -2.0 * kweight * ct.vel * sp->kbnorm[ip] *
stptr->occupation * forcesR[3];
iptr->force[ct.fpt[0]][0] -= t1 * forcesR[0];
iptr->force[ct.fpt[0]][1] -= t1 * forcesR[1];
iptr->force[ct.fpt[0]][2] -= t1 * forcesR[2];
t1 = -2.0 * kweight * ct.vel * sp->kbnorm[ip] *
stptr->occupation * forcesI[3];
iptr->force[ct.fpt[0]][0] -= t1 * forcesI[0];
iptr->force[ct.fpt[0]][1] -= t1 * forcesI[1];
iptr->force[ct.fpt[0]][2] -= t1 * forcesI[2];
#endif
} /* end for */
#if !GAMMA_PT
release_mem(tmp_psiI);
#endif
release_mem(rx);
} /* end nlforce_s */