TABLE OF CONTENTS

1. QMD-MGDFT/quench.c

QMD-MGDFT/quench.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 quench(STATE *states, REAL *vxc, REAL *vh, REAL *vnuc, 
               REAL *rho, REAL *rhocore, REAL *rhoc)
   For a fixed atomic configuration, quench the electrons to find 
   the minimum total energy 

   states: point to orbital structure (see md.h)
   vxc:    exchange correlation potential
   vh:     Hartree potential
   vnuc:   Pseudopotential 
   rho:    total valence charge density
   rhocore: core chare density only for non-linear core correction
   rhoc:   Gaussian compensating charge density

   states, vxc, vh, rho are updated

   cdfastrlx.c fastrlx.c md.c

   scf.c force.c get_te.c subdiag.c get_ke.c

     
    #include <float.h>
    #include <math.h>
    #include <stdlib.h>
    #include <stdio.h>
    #include "md.h"
    
    
    void quench(STATE *states, REAL *vxc, REAL *vh, REAL *vnuc, REAL *rho, 
                               REAL *rhocore, REAL *rhoc)
    {
    
      int outcount = 0, kpt;
      static int CONVERGENCE = FALSE;
      int st; 
      REAL KE;
    
      for(ct.steps = 0;ct.steps < ct.scfpermd;ct.steps++) {
    
        /* Perform a single self-consistent step */ 
        if (! CONVERGENCE) 
           scf(states, vxc, vh, vnuc, rho, rhocore, rhoc, &CONVERGENCE);
    
        if( !ct.steps )CONVERGENCE = FALSE;
    
    
        if (CONVERGENCE) {
    
          if(pct.thispe == 0)
            printf ("\n\n     convergence has been achieved. stopping ...\n");
    
          /* Compute the forces */
          force(rho, rhoc, vh, vxc, states);
          /* Output the forces */
          if(pct.thispe == 0)write_force();
          break;
    
        }/* end if */
    
    
        /* Check if we need to output intermediate results */
        if(outcount >= ct.outcount) {
    
    
          /* Get the total energy */
          get_te(rho, rhocore, rhoc, vh, vxc, states);
          if(pct.thispe==0)
                  printf("\n TOTAL ENERGY = %14.7f\n", ct.TOTAL); 
          
          
          /* Output the energies */
          if(pct.thispe == 0) output(states, ct.steps);
          
          
          outcount = 0;
          
    
        } /* end if */
    
    
        /* Do diagonalizations if requested */
        if(ct.diag) {
          
          if(((ct.steps % ct.diag) == 0) && (ct.end_diag>ct.steps)) {
    
            if(ct.steps) {
                    
              for(kpt=0;kpt < ct.num_kpts;kpt++)
                      subdiag(ct.kp[kpt].kstate, vh, vnuc, vxc); 
    
            } /* end if */
    
          } /* end if */
    
        } /* end if */
    
        if(pct.thispe == 0)write_occ(states);
        
        outcount++;
        
      } /* end for */
    
    
      /* Compute kinetic energy of all states */
      for(st = 0;st < ct.num_states;st++) {
        KE = get_ke(&states[st], 0);
       if(pct.thispe ==0) printf("Kinetic energy for state %d = %14.6f\n", st, KE);
      }
     
    } /* end quench */