TABLE OF CONTENTS

1. QMD-MGDFT/nlforce.c

QMD-MGDFT/nlforce.c

   Ab initio real space code with multigrid acceleration
   Quantum molecular dynamics package.
   Version: 2.1.5

   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

   void nlforce(STATE *states, int tid)
   Calculates contribution to the ionic forces from the non-local
   potentials.

   states: point to the orbital structure (see md.h_
   tid:   thread ID

   forces are added in the structure ct.ion

   force.c

   get_index.c nlforce_s.c nlforce_p.c nlforce_d.c 

    #include <float.h>
    #include <math.h>
    #include <stdio.h>
    #include <stdlib.h>
    #include "md.h"
    
    
    
    void nlforce(STATE *states, int tid)
    {
    
        int kpt, ion, idx;
        int ix, iy, iz, ip;
        int Aix[NX_GRID], Aiy[NY_GRID], Aiz[NZ_GRID];
        int *pvec, *dvec;
        int ilow, jlow, klow, ihi, jhi, khi, map, icount;
        int alloc, icut, itmp, icenter;
        SPECIES *sp;
        ION *iptr;
        int size, sizr, nion, ion_index[MAX_IONS];
    
    
    
    
        /* Grab some memory for temporary storage */
        alloc = ct.max_nlpoints;
        if(alloc < P0_BASIS)alloc = P0_BASIS;
        pvec = (int *)get_mem( 2 * alloc );
        dvec = pvec + alloc;
    
    
    
    
        /* Loop over ions */
    #if ((MACHINE_TYPE == PARALLEL_SMP) || (MACHINE_TYPE == SERIAL))
        size = ct.num_ions / ct.num_threads;
        sizr = ct.num_ions % ct.num_threads;
        for(idx = 0;idx < size;idx++) {
    
          ion_index[idx] = tid*size + idx;
    
        } /* end for */
    
        if(tid < sizr) {
          ion_index[size] = ct.num_threads*size + tid;
          size++;
        } /* end if */
    
        for(nion = 0;nion < size;nion++) {
          ion = ion_index[nion];
    
    #else
    
        for(ion = 0;ion < ct.num_ions;ion++) {
    
    #endif
    
    
            for(idx = 0;idx < alloc;idx++) {
              pvec[idx] = dvec[idx] = 0;
            }
    
            /* Generate ion pointer */
            iptr = &ct.ions[ion];
    
    
            /* Get species type */
            sp = ct.sp[iptr->species];
    
            icenter = sp->nldim / 2;
            icut = (icenter + 1) * (icenter + 1);
    
            /* Determine mapping indices or even if a mapping exists */
            map = get_index(iptr, Aix, Aiy, Aiz, &ilow, &ihi, &jlow, &jhi, &klow, &khi,
                      sp->nldim, PX0_GRID, PY0_GRID, PZ0_GRID,
                      ct.psi_nxgrid, ct.psi_nygrid, ct.psi_nzgrid,
                      &iptr->nlxcstart, &iptr->nlycstart, &iptr->nlzcstart);
    
    
            my_barrier();
    
    
            icount = idx = 0;
            /* Generate index arrays */
            for(ix = 0;ix < sp->nldim;ix++) {
    
                for(iy = 0;iy < sp->nldim;iy++) {
    
                    for(iz = 0;iz < sp->nldim;iz++) {
    
                        dvec[idx] = FALSE;
                        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[icount] =
                                   PY0_GRID * PZ0_GRID * (Aix[ix] % PX0_GRID) +
                                   PZ0_GRID * (Aiy[iy] % PY0_GRID) +
                                   (Aiz[iz] % PZ0_GRID);
    
                              dvec[idx] = TRUE;
           
                              icount++;
    
                            } /* end if */
    
                        } /* end if */
                        idx++;
    
                    } /* end for */
    
                } /* end for */
    
            } /* end for */
    
    
            /* Loop over k-points */
            for(kpt = 0;kpt < ct.num_kpts;kpt++) {
    
                /* Loop over radial projectors */
                for(ip = 0;ip < sp->num_potentials;ip++) {
    
    
                    /* Skip the local potential */
                    if(sp->lval[ip] != sp->local) {
    
                        switch(sp->lval[ip]) {
    
                            case S_STATE:
    
                                nlforce_s(ion, ip, icount, ct.kp[kpt].kstate, dvec, pvec);
                                break;
    
                            case P_STATE:
    
                                nlforce_p(ion, ip, icount, ct.kp[kpt].kstate, dvec, pvec);
                                break;
    
                            case D_STATE:
                                nlforce_d(ion, ip, icount, ct.kp[kpt].kstate, dvec, pvec);
                                break;
    
                            case F_STATE:
    
                            default:
                                error_handler("nlforce",
                                          "Angular momentum state not programmed");
    
                        } /* end switch */
     
                    } /* end if */
    
                } /* end for */
    
                my_barrier();
    
            } /* end loop over k-points */
    
    #if ((MACHINE_TYPE == PARALLEL_SMP) || (MACHINE_TYPE == SERIAL)) 
        } /* end loop over ions */
    #else
        } /* end loop over ions */
    #endif
    
        /* Release temporary memory */
        release_mem(pvec);
    
    
    } /* end nlforce */