TABLE OF CONTENTS

1. QMD-MGDFT/init_wflcao.c

QMD-MGDFT/init_wflcao.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 init_wflcao(STATE *states)
   Generates initial wavefunctions from pseudoatomic orbitals.
   if ct.occ_flag =/ 0, start with each state equally occupied

   nothing

   wave function are initialized and store in structure STATE 

   init.c

   pe2xyz.c scatter_psi.c to_cartesian.c norm_psi.c

    #include <float.h>
    #include <math.h>
    #include <stdio.h>
    #include <stdlib.h>
    #include "md.h"
    #include "inline.h"
    
    
    
    void init_wflcao(STATE *states)
    {
    
        int ion, idx, istate, sbasis;
        int ix, iy, iz, kpt;
        int pex, pey, pez;
        REAL r, ax[3], bx[3], xc, yc, zc, t1, t2, s0, px, py, pz, invdr;
        REAL *rand_ion, kdr;
        REAL *tmp_psiR, *tmp_psiI, coskdr, sinkdr;
        REAL c0, c1;
        long idum;
    
        SPECIES *sp;
        STATE *st;
        ION *iptr;
    
        c0 = sqrt(1.0 / (4.0 * PI));
        c1 = sqrt(3.0 / (4.0 * PI));
        c0 = 1.0 / c0;
        c1 = 1.0 / c1;
    c0 = 1.0;
    c1 = 1.0;
        kpt = states->kidx;
        pe2xyz(pct.thispe, &pex, &pey, &pez);
    
        /* Initialize the random number generator */
        idum = 3356;
        rand0(&idum);
    
        sbasis = (PX0_GRID + 2) * (PY0_GRID + 2) * (PZ0_GRID + 2);
        tmp_psiR = get_mem(sbasis);
        tmp_psiI = get_mem(sbasis);
        rand_ion = get_mem(4*MAX_STATES*MAX_IONS);
    
             
    
        /* Loop over states and zero out the wavefunctions */
        for(idx = 0;idx < P0_BASIS;idx++) tmp_psiR[idx] = 0.0;
        for(istate = 0;istate < ct.num_states;istate++) {
    
          st = &states[istate];
          scatter_psi(tmp_psiR, NULL, st, 0);
    #if !GAMMA_PT
          scatter_psi(NULL, tmp_psiR, st, 0);
    #endif
          if (ct.occ_flag && (ct.runflag == 0)){
            st->occupation = ct.nel/ct.num_states;
          }
    
        } /* end for */
    
        for(idx = 0;idx < ct.num_states*ct.num_ions*4;idx++) {
          rand_ion[idx] = 1.0 + 0.2*(rand0(&idum) - 0.5);
        } /* 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];
            invdr = ONE / sp->drnlig;
    
    
            /* Loop over the grid */
            idx = 0;
            xc = ((REAL)pex) / PE_X;
            for(ix = 0;ix < PX0_GRID;ix++) {
    
                yc = ((REAL)pey) / PE_Y;
                for(iy = 0;iy < PY0_GRID;iy++) {
    
                    zc = ((REAL)pez) / PE_Z;
                    for(iz = 0;iz < PZ0_GRID;iz++) {
    
                      ax[0] = xc - iptr->xtal[0];
                      ax[1] = yc - iptr->xtal[1];
                      ax[2] = zc - iptr->xtal[2];
    
                      r = metric(ax) + 1.0e-20;
    
    
                      /* s-orbital */
                      s0 = linint(&sp->psilig[0][0], r, invdr);
    
                      /* p-orbital */
                      t1 = linint(&sp->psilig[1][0], r, invdr);
    
                      to_cartesian(ax, bx);
    
    #if !GAMMA_PT
                      kdr = ct.kp[kpt].kvec[0] * bx[0] +
                      ct.kp[kpt].kvec[1] * bx[1] +
                      ct.kp[kpt].kvec[2] * bx[2];                           
    
                      coskdr = cos(kdr);
                      sinkdr = sin(kdr);
    #endif
    
                      px = bx[0] * t1 / r;
                      py = bx[1] * t1 / r;
                      pz = bx[2] * t1 / r;
    
                      /* Loop over states */
                      st = &states[0];
    #if GAMMA_PT
                      st->psiR[idx] += s0;  /* S-component only in lowest state */
    #else
                      st->psiR[idx] += s0; 
                      st->psiI[idx] += 0.0;
    #endif
                      for(istate = 1;istate < ct.num_states;istate++) {
    
                          st = &states[istate];
    #if GAMMA_PT
                          st->psiR[idx] += c0*rand_ion[ion*istate]*s0 + 
                                         c1*rand_ion[ion*istate+1]*px + 
                                         c1*rand_ion[ion*istate+2]*py + 
                                         c1*rand_ion[ion*istate+3]*pz;
    #else
    
                          t2 = c0*rand_ion[ion*istate]*s0 +
                               c1*rand_ion[ion*istate+1]*px +
                               c1*rand_ion[ion*istate+2]*py +
                               c1*rand_ion[ion*istate+3]*py;
                          st->psiR[idx] += coskdr * t2;
                          st->psiI[idx] += sinkdr * t2;
    #endif
    
                      } /* end for */
    
    
                      idx++;
                      zc += ct.hzgrid;
                    } /* end for */
    
                    yc += ct.hygrid;
                } /* end for */
    
                xc += ct.hxgrid;
            } /* end for */
    
        } /* end for */
    
    
        /* Normalize the orbitals */
        for(istate = 0;istate < ct.num_states;istate++) {
          st = &states[istate];
          norm_psi(st, ONE);
        } /* end for */
    
        release_mem(rand_ion);
        release_mem(tmp_psiI);
        release_mem(tmp_psiR);
    
    } /* end init_wflcao */