TABLE OF CONTENTS

1. QMD-MGDFT/fastrlx.c

QMD-MGDFT/fastrlx.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 fastrlx(STATE *states, REAL *vxc, REAL *vh, REAL *vnuc,
                REAL *rho, REAL *rhocore, REAL *rhoc)
   drive routine for fast relax.

   states: all wave functions (see md.h)
   vxc: exchange-correlation potential
   vh:  Hartree potential
   vnuc: pseudopotential
   rho:  total charge density
   rhocore: charge density of core electrons, only useful when we 
            include non-linear core correction for pseudopotential.
   rhoc:    compensating charge density

   all the above inputs are updated

   md.c

   quench.c to_crystal.c to_cartesian.c init_nuc.c get_nlop.c scf.c
   get_te.c force.c md_fastrelax.c

    #include <float.h>
    #include <math.h>
    #include <stdlib.h>
    #include <stdio.h>
    #include "md.h"
    
    
    
    /* Local function prototypes */
    void md_fastrelax(void);
    void movie(FILE *);
    
    
    
    void fastrlx(STATE *states, REAL *vxc, REAL *vh, REAL *vnuc,
             REAL *rho, REAL *rhocore, REAL *rhoc)
    {
    
      int steps, it, ion;
      int CONVERGENCE = FALSE;
      int CONV_FORCE = FALSE;
      int temp_num;
      int diagcount = 0;
      ION *iptr;
      FILE *mfp=NULL;
    
      if (ct.override_atoms==1) {
        quench(states, vxc, vh, vnuc, rho, rhocore, rhoc);
      }
    
      /* open movie file */
      if(pct.thispe == 0) {
        mfp = fopen("traj.rmv","w");
        if(setvbuf(mfp,(char *)NULL,_IOFBF,4096*16) != 0)
          printf("\n Warning: cant allocate movie io buffer size\n");
    
        /* output initial frame */
        movie(mfp);
      }
    
    
      for(ct.steps = 0;ct.steps < ct.num_steps;ct.steps++) {
    
          /* enforce periodic boundary conditions on the ions */
          for(it = 0;it < ct.num_ions;it++){
    
            iptr = &ct.ions[it];
    
            /* to_crystal enforces periodic boundary conditions */
            to_crystal(iptr->xtal,iptr->crds);
            to_cartesian(iptr->xtal, iptr->crds);
    
          } /* end for */
    
    
          CONVERGENCE = FALSE;
          
          for(steps = 0;steps < ct.scfpermd;steps++) {
    
              /* Perform a single self-consistent step */
            if (! CONVERGENCE) {
              
              scf(states, vxc, vh, vnuc, rho, rhocore, rhoc, &CONVERGENCE);
             
    
              /* Do diagonalizations if requested */
              if(ct.diag) {
                if(diagcount > ct.diag) {
                  subdiag(states, vh, vnuc, vxc);
                  diagcount = 0;
                } /* end if */
                diagcount++;
              } /* end if */
              
              /* Get the total energy */
              get_te(rho, rhocore, rhoc, vh, vxc, states);
              
            } /* end if "convergence" */
            
          } /* end for steps */ 
    
    
          /* Compute the forces */
          force(rho, rhoc, vh, vxc, states);
    
          /* Output the forces */
          if(pct.thispe == 0)write_force();
    
          /* out put the milliken charges */
          if(ct.domilliken == 1) {
            ct.milliken = 1;
            get_milliken(states);
            get_nlop();
          } /* if(ct.domilliken == 1) */
    
          CONV_FORCE = TRUE;
          for(ion = 0;ion < ct.num_ions;ion++) {
            if(ct.ions[ion].movable){
              CONV_FORCE = (CONV_FORCE &&
                            ( (fabs(ct.ions[ion].force[ct.fpt[0]][0]) 
                               < ct.thr_frc) &&
                              (fabs(ct.ions[ion].force[ct.fpt[0]][1]) 
                               < ct.thr_frc) &&
                              (fabs(ct.ions[ion].force[ct.fpt[0]][2]) 
                               < ct.thr_frc)));
            } /* end if "movable ions" */
          } /* end for ion */
    
          if (! CONV_FORCE) {
            
            /* Move the ions */
            md_fastrelax();
    
            /* Update items that change when the ionic coordinates change */
            init_nuc(vnuc, rhoc, rhocore);
            get_nlop();
    
            /* write out frame */
            if(pct.thispe == 0) {
              movie(mfp);
            }
          }
          else{
            /* close moviefile */
            if(pct.thispe == 0) {
              printf ("\n\n     convergence in forces has been achieved. stopping ...\n");
    
              fclose(mfp);
            }
            
            return;
          }
            
            
      } /* end for */
    
      /* close moviefile */
      if(pct.thispe == 0) {
        fclose(mfp);
      }
     
    } /* end fastrlx */