Generated from nlforce_d.c with ROBODoc v3.2.2 on Fri Sep 14 14:23:20 2001
TABLE OF CONTENTS
- QMD-MGDFT/nlforce_d.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 d-projectors.
INPUTS
ion: ion index
ip: momentum l = 2
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
/*
nlforce_p.c
*/
#include <float.h>
#include <math.h>
#include <stdio.h>
#include <stdlib.h>
#include "md.h"
#include "inline.h"
/* Determines forces from a set of d-projectors */
void nlforce_d(int ion, int ip, int icount, STATE *states, int *dvec, int *pvec)
{
#if 0
int idx, ix, iy, iz, docount, st, incx = 1;
REAL r, rsq, rsqd, r3, x, y, z, xc, yc, zc;
REAL t1, t2, t3, invdr;
REAL fx, fy, fz;
REAL *tmp_psi, *loc_psi;
REAL *prjptr_1, *prjptr_2, *prjptr_3, *prjptr_4, *prjptr_5;
REAL *r1_x, *r1_y, *r1_z;
REAL *r2_x, *r2_y, *r2_z;
REAL *r3_x, *r3_y, *r3_z;
REAL *r4_x, *r4_y, *r4_z;
REAL *r5_x, *r5_y, *r5_z;
REAL dy1_dx, dy1_dy, dy1_dz;
REAL dy2_dx, dy2_dy, dy2_dz;
REAL dy3_dx, dy3_dy, dy3_dz;
REAL dy4_dx, dy4_dy, dy4_dz;
REAL dy5_dx, dy5_dy, dy5_dz;
REAL dt2_dx, dt2_dy, dt2_dz;
REAL ax[3], bx[3], y1, y2, y3, y4, y5;
REAL forces[4];
SPECIES *sp;
ION *iptr;
STATE *stptr;
r1_x = get_mem(25 * P0_BASIS);
r1_y = r1_x + P0_BASIS;
r1_z = r1_y + P0_BASIS;
r2_x = r1_z + P0_BASIS;
r2_y = r2_x + P0_BASIS;
r2_z = r2_y + P0_BASIS;
r3_x = r2_z + P0_BASIS;
r3_y = r3_x + P0_BASIS;
r3_z = r3_y + P0_BASIS;
r4_x = r3_z + P0_BASIS;
r4_y = r4_x + P0_BASIS;
r4_z = r4_y + P0_BASIS;
r5_x = r4_z + P0_BASIS;
r5_y = r5_x + P0_BASIS;
r5_z = r5_y + P0_BASIS;
tmp_psi = r5_z + P0_BASIS;
loc_psi = tmp_psi + P0_BASIS;
prjptr_1 = loc_psi + P0_BASIS;
prjptr_2 = prjptr_1+ P0_BASIS;
prjptr_3 = prjptr_2+ P0_BASIS;
prjptr_4 = prjptr_3+ P0_BASIS;
prjptr_5 = prjptr_4+ P0_BASIS;
/* Get ion pointer */
iptr = &ct.ions[ion];
/* Get species type */
sp = ct.sp[iptr->species];
invdr = ONE / sp->drnlig;
t3 = sqrt(3.0);
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);
x = bx[0];
y = bx[1];
z = bx[2];
rsq = x*x + y*y + z*z;
r = sqrt(rsq);
r3 = r * rsq + 1.0e-10;
rsqd = rsq + 1.0e-10;
t1 = t3 * linint(&sp->drkbplig[ip][0], r, invdr);
t2 = t3 * linint(&sp->kbplig[ip][0], r, invdr);
r += 1.0e-10;
y1 = x * y / rsqd;
y2 = y * z / rsqd;
y3 = z * x / rsqd;
y4 = (x*x - y*y) / (2.0*rsqd);
y5 = (t3*z*z - rsq/t3) / (2.0*rsqd);
dy1_dx = (y - 2.0*x*y/rsqd) / rsqd;
dy1_dy = (x - 2.0*x*y/rsqd) / rsqd;
dy1_dz = -2.0*x*y*z/r3;
dy2_dx = -2.0*x*y*z/r3;
dy2_dy = (z - 2.0*z*y/rsqd) / rsqd;
dy2_dz = (y - 2.0*z*y/rsqd) / rsqd;
dy3_dx = (z - 2.0*z*x/rsqd) / rsqd;
dy3_dy = -2.0*x*y*z/r3;
dy3_dz = (x - 2.0*z*x/rsqd) / rsqd;
dy4_dx = x*(1 - x*(x*x - y*y)/rsqd)/rsqd;
dy4_dy = y*(-1 - y*(x*x - y*y)/rsqd)/rsqd;
dy4_dz = -z*(x*x - y*y) / (rsqd * rsqd);
dy5_dx = - t3*x*z*z / (rsqd * rsqd);
dy5_dy = - t3*y*z*z / (rsqd * rsqd);
dy5_dz = t3 * (z - z*z/rsqd) / rsqd;
dt2_dx = x * t1 / r;
dt2_dy = y * t1 / r;
dt2_dz = z * t1 / r;
prjptr_1[docount] = t2 * y1;
prjptr_2[docount] = t2 * y2;
prjptr_3[docount] = t2 * y3;
prjptr_4[docount] = t2 * y4;
prjptr_5[docount] = t2 * y5;
r1_x[docount] = t2 * dy1_dx + y1 * dt2_dx;
r1_y[docount] = t2 * dy1_dy + y1 * dt2_dy;
r1_z[docount] = t2 * dy1_dz + y1 * dt2_dz;
r2_x[docount] = t2 * dy2_dx + y2 * dt2_dx;
r2_y[docount] = t2 * dy2_dy + y2 * dt2_dy;
r2_z[docount] = t2 * dy2_dz + y2 * dt2_dz;
r3_x[docount] = t2 * dy3_dx + y3 * dt2_dx;
r3_y[docount] = t2 * dy3_dy + y3 * dt2_dy;
r3_z[docount] = t2 * dy3_dz + y3 * dt2_dz;
r4_x[docount] = t2 * dy4_dx + y4 * dt2_dx;
r4_y[docount] = t2 * dy4_dy + y4 * dt2_dy;
r4_z[docount] = t2 * dy4_dz + y4 * dt2_dz;
r5_x[docount] = t2 * dy5_dx + y5 * dt2_dx;
r5_y[docount] = t2 * dy5_dy + y5 * dt2_dy;
r5_z[docount] = t2 * dy5_dz + y5 * dt2_dz;
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 */
gather_psi(tmp_psi, stptr, 0);
for(idx = 0;idx < icount;idx++) loc_psi[idx] = tmp_psi[pvec[idx]];
forces[0] = ct.vel * sdot(&icount, r1_x, &incx, loc_psi, &incx);
forces[1] = ct.vel * sdot(&icount, r1_y, &incx, loc_psi, &incx);
forces[2] = ct.vel * sdot(&icount, r1_z, &incx, loc_psi, &incx);
forces[3] = ct.vel * sdot(&icount, prjptr_1, &incx, loc_psi, &incx);
idx = 4;
global_sums(forces, &idx);
t1 = -2.0 * forces[3] * iptr->pd[ip] * stptr->occupation;
if(pct.thispe == 0)
printf("a %d %d %f %f\n", ion, st, t1, forces[0]);
fflush(NULL);
iptr->force[ct.fpt[0]][0] -= t1 * forces[0];
iptr->force[ct.fpt[0]][1] -= t1 * forces[1];
iptr->force[ct.fpt[0]][2] -= t1 * forces[2];
forces[0] = ct.vel * sdot(&icount, r2_x, &incx, loc_psi, &incx);
forces[1] = ct.vel * sdot(&icount, r2_y, &incx, loc_psi, &incx);
forces[2] = ct.vel * sdot(&icount, r2_z, &incx, loc_psi, &incx);
forces[3] = ct.vel * sdot(&icount, prjptr_2, &incx, loc_psi, &incx);
idx = 4;
global_sums(forces, &idx);
t1 = -2.0 * forces[3] * iptr->pd[ip] * stptr->occupation;
if(pct.thispe == 0)
printf("b %d %d %f %f\n", ion, st, t1, forces[0]);
fflush(NULL);
iptr->force[ct.fpt[0]][0] -= t1 * forces[0];
iptr->force[ct.fpt[0]][1] -= t1 * forces[1];
iptr->force[ct.fpt[0]][2] -= t1 * forces[2];
forces[0] = ct.vel * sdot(&icount, r3_x, &incx, loc_psi, &incx);
forces[1] = ct.vel * sdot(&icount, r3_y, &incx, loc_psi, &incx);
forces[2] = ct.vel * sdot(&icount, r3_z, &incx, loc_psi, &incx);
forces[3] = ct.vel * sdot(&icount, prjptr_3, &incx, loc_psi, &incx);
idx = 4;
global_sums(forces, &idx);
t1 = -2.0 * forces[3] * iptr->pd[ip] * stptr->occupation;
if(pct.thispe == 0)
printf("c %d %d %f %f\n", ion, st, t1, forces[0]);
fflush(NULL);
iptr->force[ct.fpt[0]][0] -= t1 * forces[0];
iptr->force[ct.fpt[0]][1] -= t1 * forces[1];
iptr->force[ct.fpt[0]][2] -= t1 * forces[2];
forces[0] = ct.vel * sdot(&icount, r4_x, &incx, loc_psi, &incx);
forces[1] = ct.vel * sdot(&icount, r4_y, &incx, loc_psi, &incx);
forces[2] = ct.vel * sdot(&icount, r4_z, &incx, loc_psi, &incx);
forces[3] = ct.vel * sdot(&icount, prjptr_4, &incx, loc_psi, &incx);
idx = 4;
global_sums(forces, &idx);
t1 = -2.0 * forces[3] * iptr->pd[ip] * stptr->occupation;
if(pct.thispe == 0)
printf("d %d %d %f %f\n", ion, st, t1, forces[0]);
fflush(NULL);
iptr->force[ct.fpt[0]][0] -= t1 * forces[0];
iptr->force[ct.fpt[0]][1] -= t1 * forces[1];
iptr->force[ct.fpt[0]][2] -= t1 * forces[2];
forces[0] = ct.vel * sdot(&icount, r5_x, &incx, loc_psi, &incx);
forces[1] = ct.vel * sdot(&icount, r5_y, &incx, loc_psi, &incx);
forces[2] = ct.vel * sdot(&icount, r5_z, &incx, loc_psi, &incx);
forces[3] = ct.vel * sdot(&icount, prjptr_5, &incx, loc_psi, &incx);
forces[3] = ct.vel * sdot(&icount, prjptr_4, &incx, loc_psi, &incx);
idx = 4;
global_sums(forces, &idx);
t1 = -2.0 * forces[3] * iptr->pd[ip] * stptr->occupation;
if(pct.thispe == 0)
printf("e %d %d %f %f\n", ion, st, t1, forces[0]);
fflush(NULL);
iptr->force[ct.fpt[0]][0] -= t1 * forces[0];
iptr->force[ct.fpt[0]][1] -= t1 * forces[1];
iptr->force[ct.fpt[0]][2] -= t1 * forces[2];
} /* end for */
release_mem(r1_x);
#endif
} /* end nlforce_d */