31. Pseudopotential Plane-Wave DFT (PSPW)

A pseudopotential plane-wave (PSPW) module, which can perform Car-Parrinello simulations, is being implemented into the NWChem program package. This module complements the capabilities of the more traditional Gaussian function based approaches by including code which allows for the calculation of density functional theory total energies and forces with the technology based on plane-wave basis sets and pseudopotentials. Consistent with NWChem's philosophy this module is able to run on a variety of architectures, including parallel supercomputers.

The advantage of a PSPW method is that it has been shown to have an accuracy close to chemical accuracy for many applications, yet is still fast enough to treat systems containing hundreds of atoms. Another significant advantage is its ability to simulate dynamics on a ground state potential surface directly at run-time. This method's efficiency and accuracy make it a desirable first principles method of simulation in the study of complex molecular, liquid, and solid state systems. Applications for this first principles method include the calculation of free energies, search for global minima, explicit simulation of solvated molecules, and simulations of complex vibrational modes that cannot be described within the harmonic approximation.

Section 31.1 describes the tasks contained within the PSPW module. The RTDB entries and datafiles used by the PSPW module are described in section 31.2. Car-Parrinello output data files are described in section 31.2.9, and the minimization and Car-Parrinello algorithms are described in sections 31.3-31.4. Examples of how to setup and run a steepest descent simulation and a Car-Parrinello simulation are presented in sections 31.5-31.6. Finally in section 31.7 the capabilities and limitations of the PSPW module are discussed.

If you are a first time user of this module it is recommended that you skip the next five sections and proceed directly to the tutorials in sections 31.5-31.6.

31.1 PSPW Tasks

All input to the PSPW Tasks is contained within the compound PSPW block,

PSPW
   ...
END

To perform an actual calculation a TASK PSPW directive is used. The TASK cg_pspw directive may also be used (see below). The format for the TASK PSPW directive is the TASK directive followed by the PSPW string, and after that an $<$operation$>$ string is required. The TASK PSPW directive for PSPW calculations is of the following form:

TASK PSPW [steepest_descent           ||
           conjugate_gradient         ||
           Car-Parrinello             ||
           psp_formatter              ||
           wavefunction_initializer   ||
           v_wavefunction_initializer ||
           wavefunction_expander      ||
           psp_generator              ||
           (no default)]

Currently available tasks are listed. Note that unlike most NWChem modules, the PSPW module does not contain an energy operation. This means that there is no default operation and hence an operation must always be specified.

PSPW tasks are of two types. The first type of task is used to setup or change needed data files for a PSPW simulation. Tasks of this type are psp_formatter, wavefunction_initializer, v_wavefunction_initializer, wavefunction_expander, and psp_generator. The second type of task is used to actually run a PSPW simulation. Tasks of this type are steepest_descent, and Car-Parrinello. The following subsections describe the input to these tasks.

The PSPW module is also interfaced to driver, stepper, and freqency task operations through the cg_pspw theory directive, see section 5.10.1. The PSPW conjugate_gradient sub-module is used in this interface. Specifically the following task operations are interafaced using the cg_pspw theory directive.

TASK cg_pspw energy      #equivalent to TASK PSPW conjugate_gradient           
TASK cg_pspw gradient         
TASK cg_pspw optimize         
TASK cg_pspw saddle           
TASK cg_pspw freqencies       
TASK cg_pspw vib

31.1.1 STEEPEST_DESCENT

The steepest_descent task is used to optimize the one-electron orbitals with respect to the total energy. In addition it can also be used to optimize geometries. This method is meant to be used for coarse optimization of the one-electron orbitals. A detailed description of the this method is described in section 31.3.1

Input to the steepest_descent simulation is contained within the steepest_descent sub-block.

PSPW
  ...
  STEEPEST_DESCENT
     ...
  END
  ...
END

To run a steepest_descent calculation the following directive is used:

TASK PSPW steepest_descent

The steepest_descent sub-block contains a great deal of input, including pointers to data, as well as parameter input. Listed below is the format of a STEEPEST_DESCENT sub-block.

PSPW
...
   STEEPEST_DESCENT
      CELL_NAME: <string cell_name>
      [GEOMETRY_OPTIMIZE]
      (FORMATTED_FILENAME: <string formatted_name> \
       ...)
      INPUT_WAVEFUNCTION_FILENAME:  <string input_wavefunctions  default input_movecs>
      OUTPUT_WAVEFUNCTION_FILENAME: <string output_wavefunctions default input_movecs>
      FAKE_MASS: <real fake_mass default 400000.0>
      TIME_STEP: <real time_step default 5.8>
      LOOP: <integer inner_iteration outer_iteration default 10 1>
      TOLERANCES: <real tole tolc tolr default 1.0d-9 1.0d-9 1.0d-4>
      ENERGY_CUTOFF:       <real ecut default (see input desciption)>
      WAVEFUNCTION_CUTOFF: <real wcut default (see input description)>
      EWALD_NCUT: <integer ncut default 1>
      EWALD_RCUT: <real rcut default (see input description)>
      EXCHANGE_CORRELATION: (Vosko || PBE96  default Vosko)
      [MULLIKEN]

   END
...

END

The following list describes the input for the STEEPEST_DESCENT sub-block.

• $<$cell_name$>$ - name that points to the the simulation_cell named $<$cell_name$>$. See section 31.2.3.
• GEOMETRY_OPTIMIZE - optional keyword which if specified turns on geometry optimization.
• $<$formatted_name$>$ - name that points to a formatted_pseudopotential file. A file must be specified for each kind of ion in the simulation.
• $<$input_wavefunctions$>$ - name that points to a file containing one-electron orbitals
• $<$output_wavefunctions$>$ - name that will point to file containing the one-electron orbitals at the end of the run.
• $<$fake_mass$>$ - value for the electronic fake mass ($\mu$).
• $<$time_step$>$ - value for the time step ($\Delta t$).
• $<$inner_iteration$>$ - number of iterations between the printing out of energies and tolerances
• $<$outer_iteration$>$ - number of outer iterations
• $<$tole$>$ - value for the energy tolerance.
• $<$tolc$>$ - value for the one-electron orbital tolerance.
• $<$tolr$>$ - value for the ion position tolerance.
• $<$ecut$>$ - value for the cutoff energy used to define the density. Default is set to be the maximum value that will fit within the simulation_cell $<$cell_name$>$.
• $<$wcut$>$ - value for the cutoff energy used to define the one-electron orbitals. Default is set to be the maximum value that will fit within the simulation_cell $<$cell_name$>$.
• $<$ncut$>$ - value for the number of unit cells to sum over (in each direction) for the real space part of the Ewald summation. Note Ewald summation is only used if the simulation_cell is periodic.
• $<$rcut$>$ - value for the cutoff radius used in the Ewald summation. Note Ewald summation is only used if the simulation_cell is periodic.
  Default set to be $\frac{MIN(\left\vert \vec{a_i} \right\vert)}{\pi}, i=1,2,3$.
• (Vosko $\vert\vert$ PBE96) - Choose between Vosko et al's LDA parameterization or the Perdew, Burke, and Erzherhoff GGA functional.
• MULLIKEN - optional keyword which if specified causes a Mulliken anaysis to be performed at the end of the simulation. For this option to work angular momentum weights for each kind of atom have to be entered into the RTDB using the PSPW ANALYSIS sub-block (see section 31.2.4.

31.1.2 CONJUGATE_GRADIENT

The conjugate_gradient task is used to optimize the one-electron orbitals with respect to the total energy. This method should be used for finer optimization. It is better to use the steepest_descent task when intitially optimizing the one-electron orbitals. A detailed description of the this method is described in section 31.3.2.

Input to the conjugate_gradient simulation is contained within the conjugate_gradient sub-block.

PSPW
  ...
  CONJUGATE_GRADIENT
     ...
  END
  ...
END

To run a conjugate_gradient calculation the following directive is used:

TASK PSPW conjugate_gradient

The conjugate_gradient sub-block contains a great deal of input, including pointers to data, as well as parameter input. Listed below is the format of a conjugate_gradient sub-block.

PSPW
...
   CONJUGATE_GRADIENT
      CELL_NAME: <string cell_name>
      (FORMATTED_FILENAME: <string formatted_name> \
       ...)
      INPUT_WAVEFUNCTION_FILENAME:  <string input_wavefunctions  default input_movecs>
      OUTPUT_WAVEFUNCTION_FILENAME: <string output_wavefunctions default input_movecs>
      [FAKE_MASS: <real fake_mass default 400000.0>]
      [TIME_STEP: <real time_step default 5.8>]
      LOOP: <integer inner_iteration outer_iteration default 10 1>
      TOLERANCES: <real tole tolc default 1.0e-9 1.0e-9>
      ENERGY_CUTOFF:       <real ecut default (see input description)>
      WAVEFUNCTION_CUTOFF: <real wcut default (see input description)>
      EWALD_NCUT: <integer ncut default 1>]
      EWALD_RCUT: <real rcut default (see input description)>
      EXCHANGE_CORRELATION: (Vosko || PBE96  default Vosko)
      [MULLIKEN]

   END
...

END

The following list describes the input for the CONJUGATE_GRADIENT sub-block.

• $<$cell_name$>$ - name that points to the the simulation_cell named $<$cell_name$>$. See section 31.2.3.
• $<$formatted_name$>$ - name that points to a formatted_pseudopotential file. A file must be specified for each kind of ion in the simulation.
• $<$input_wavefunctions$>$ - name that points to a file containing one-electron orbitals
• $<$output_wavefunctions$>$ - name that will point to file containing the one-electron orbitals at the end of the run.
• $<$fake_mass$>$ - value for the electronic fake mass ($\mu$). This parameter is not presently used in a conjugate gradient simulation
• $<$time_step$>$ - value for the time step ($\Delta t$). This parameter is not presently used in a conjugate gradient simulation.
• $<$inner_iteration$>$ - number of iterations between the printing out of energies and tolerances
• $<$outer_iteration$>$ - number of outer iterations
• $<$tole$>$ - value for the energy tolerance.
• $<$tolc$>$ - value for the one-electron orbital tolerance.
• $<$ecut$>$ - value for the cutoff energy used to define the density. Default is set to be the maximum value that will fit within the simulation_cell $<$cell_name$>$.
• $<$wcut$>$ - value for the cutoff energy used to define the one-electron orbitals. Default is set to be the maximum value that will fix within the simulation_cell $<$cell_name$>$.
• $<$ncut$>$ - value for the number of unit cells to sum over (in each direction) for the real space part of the Ewald summation. Note Ewald summation is only used if the simulation_cell is periodic.
• $<$rcut$>$ - value for the cutoff radius used in the Ewald summation. Note Ewald summation is only used if the simulation_cell is periodic.
  Default set to be $\frac{MIN(\left\vert \vec{a_i} \right\vert)}{\pi}, i=1,2,3$.
• (Vosko $\vert\vert$ PBE96) - Choose between Vosko et al's LDA parameterization or the Perdew, Burke, and Erzherhoff GGA functional.
• MULLIKEN - optional keyword which if specified causes a Mulliken anaysis to be performed at the end of the simulation. For this option to work angular momentum weights for each kind of atom have to be entered into the RTDB using the PSPW ANALYSIS sub-block (see section 31.2.4.

31.1.3 Car-Parrinello

The Car-Parrinello task is used to perform ab initio molecular dynamics using the scheme developed by Car and Parrinello. In this unified ab initio molecular dynamics scheme the motion of the ion cores is coupled to a fictitious motion for the Kohn-Sham orbitals of density functional theory. Constant energy or constant temperature simulations can be performed. A detailed description of this method is described in section 31.4.

Input to the Car-Parrinello simulation is contained within the Car-Parrinello sub-block.

PSPW
  ...
  Car-Parrinello
     ...
  END
  ...
END

To run a Car-Parrinello calculation the following directive is used:

TASK PSPW Car-Parrinello

The Car-Parrinello sub-block contains a great deal of input, including pointers to data, as well as parameter input. Listed below is the format of a Car-Parrinello sub-block.

PSPW
...
   Car-Parrinello
      CELL_NAME: <string cell_name>
      (FORMATTED_FILENAME: <string formatted_name> \
       ...)
      INPUT_WAVEFUNCTION_FILENAME:    <string input_wavefunctions    default input_movecs>
      OUTPUT_WAVEFUNCTION_FILENAME:   <string output_wavefunctions   default input_movecs>
      INPUT_V_WAVEFUNCTION_FILENAME:  <string input_v_wavefunctions  default input_vmovecs>
      OUTPUT_V_WAVEFUNCTION_FILENAME: <string output_v_wavefunctions default input_vmovecs>
      FAKE_MASS: <real fake_mass default default 1000.0>
      TIME_STEP: <real time_step default 5.0>
      LOOP: <integer inner_iteration outer_iteration default 10 1>
      SCALING: <real scale_c scale_r default 1.0 1.0>
      ENERGY_CUTOFF:       <real ecut default (see input description)>
      WAVEFUNCTION_CUTOFF: <real wcut default (see input description)>
      EWALD_NCUT: <integer ncut default 1>
      EWALD_RCUT: <real rcut    default (see input description)>
      EXCHANGE_CORRELATION: (Vosko || PBE96  default Vosko)
      [Nose-Hoover: <real Period_electron Temperature_electrion Period_ion Temperature_ion 
                          default 100.0 298.15 100.0 298.15>]
      XYZ_FILENAME: <string xyz_filename default XYZ>
      EMOTION_FILENAME: <string emotion_filename default EMOTION>
      HMOTION_FILENAME: <string hmotion_filename default HMOTION>
      OMOTION_FILENAME: <string omotion_filename default OMOTION>
      EIGMOTION_FILENAME: <string eigmotion_filename default EIGMOTION>
      ION_MOTION_FILENAME: <string ion_motion_filename default MOTION>

   END
...

END

The following list describes the input for the Car-Parrinello sub-block.

• $<$cell_name$>$ - name that points to the the simulation_cell named $<$cell_name$>$. See section 31.2.3.
• $<$formatted_name$>$ - name that points to a formatted_pseudopotential file. A file must be specified for each kind of ion in the simulation.
• $<$input_wavefunctions$>$ - name that points to a file containing one-electron orbitals
• $<$output_wavefunctions$>$ - name that will point to file containing the one-electron orbitals at the end of the run.
• $<$input_v_wavefunctions$>$ - name that points to a file containing one-electron orbital velocities.
• $<$output_v_wavefunctions$>$ - name that will point to file containing the one-electron orbital velocities at the end of the run.
• $<$fake_mass$>$ - value for the electronic fake mass ($\mu$).
• $<$time_step$>$ - value for the Verlet integration time step ($\Delta t$).
• $<$inner_iteration$>$ - number of iterations between the printing out of energies.
• $<$outer_iteration$>$ - number of outer iterations
• $<$scale_c$>$ - value for the initial velocity scaling of the one-electron orbital velocities.
• $<$scale_r$>$ - value for the initial velocity scaling of the ion velocities.
• $<$ecut$>$ - value for the cutoff energy used to define the density. Default is set to be the maximum value that will fit within the simulation_cell $<$cell_name$>$.
• $<$wcut$>$ - value for the cutoff energy used to define the one-electron orbitals. Default is set to be the maximum value that will fit within the simulation_cell $<$cell_name$>$.
• $<$ncut$>$ - value for the number of unit cells to sum over (in each direction) for the real space part of the Ewald summation. Note Ewald summation is only used if the simulation_cell is periodic.
• $<$rcut$>$ - value for the cutoff radius used in the Ewald summation. Note Ewald summation is only used if the simulation_cell is periodic.
  Default set to be $\frac{MIN(\left\vert \vec{a_i} \right\vert)}{\pi}, i=1,2,3$.
• (Vosko $\vert\vert$ PBE96) - Choose between Vosko et al's LDA parameterization or the Perdew, Burke, and Erzherhoff GGA functional.
• Nose-Hoover: - optional subblock which if specified causes the simulation to perform Nose-Hoover dynamics. If this subblock is not specified the simulation performs constant energy dyanmics. See section 31.4.2 for a description of the parameters.
  • $<$Period_electron$>$ $\equiv$ $P_{electron}$ - estimated period for fictitious electron thermostat.
  • $<$Temperature_electron$>$ $\equiv$ $T_{electron}$ - temperature for fictitious electron motion
  • $<$Period_ion$>$ $\equiv$ $P_{ion}$ - estimated period for ionic thermostat
  • $<$Temperature_ion$>$ $\equiv$ $T_{ion}$ - temperature for ion motion
• $<$xyz_filename$>$ - name that points to the XYZ motion file generated
• $<$emotion_filename$>$ - name that points to the emotion motion file generated. See section 31.2.9 for a description of the datafile.
• $<$hmotion_filename$>$ - name that points to the hmotion motion file generated. See section 31.2.9 for a description of the datafile.
• $<$eigmotion_filename$>$ - name that points to the eigmotion motion file generated. See section 31.2.9 for a description of the datafile.
• $<$ion_motion_filename$>$ - name that points to the ion_motion motion file generated. See section 31.2.9 for a description of the datafile.
• MULLIKEN - optional keyword which if specified causes an omotion motion file to be created. For this option to work angular momentum weights for each kind of atom have to be entered into the RTDB using the PSPW ANALYSIS sub-block (see section 31.2.4.
• $<$omotion_filename$>$ - name that points to the omotion motion file generated. See section 31.2.9 for a description of the datafile.

31.1.4 PSP_FORMATTER

The psp_formatter task takes a non-separable pseudopotential defined in one-dimension real-space in a one-dimensional psp datafile and does two things to it. First it puts it into the semi-local form suggested by Kleinman and Bylander and then it expands in in a periodic Fourier series defined by the simulation cell in the RTDB.

Input to the PSP_FORMATTER task is contained within the PSP_FORMATTER sub-block.

PSPW
  ...
  PSP_FORMATTER
     ...
  END
  ...
END

To run a PSP_FORMATTER calculation the following directive is used:

TASK PSPW PSP_FORMATTER

Listed below is the format of a PSP_FORMATTER sub-block.

PSPW
... 
   PSP_FORMATTER
      CELL_NAME:          <string cell_name> 
      PSP_FILENAME:       <string psp_name>
      FORMATTED_FILENAME: <string formatted_name>
      LOCP: (s||p||d||f||g default (see input description))
      LMAX: (s||p||d||f||g default (see input description))
   end
...
END

The following list describes the input for the PSP_FORMATTER sub-block.

• $<$cell_name$>$ - name that points to the simulation_cell named $<$cell_name$>$. See section 31.2.3.
• $<$psp_name$>$ - name that points to a one-dimensional pseudopotential datafile.
• $<$formatted_name$>$ - name that points to a formatted_pseudopotential datafile.
• LMAX (s$\vert\vert$p$\vert\vert$d$\vert\vert$f$\vert\vert$g) - used to specify the maximum number of angular potentials to use. The default is set to the maximum angular momentum in $<$psp_name$>$.
• LOCP (s$\vert\vert$p$\vert\vert$d$\vert\vert$f$\vert\vert$g) - used to specify which angular potential is to used as the local potential. The default is is set to LMAX.

31.1.5 WAVEFUNCTION_INTITIALIZER

The wavefunction_initializer task is used to generate an initial wavefunction datafile. Input to the WAVEFUNCTION_INITIALIZER task is contained within the WAVEFUNCTION_INITIALIZER sub-block.

PSPW
  ...
  WAVEFUNCTION_INITIALIZER
     ...
  END
  ...
END

To run a WAVEFUNCTION_INITIALIZER calculation the following directive is used:

TASK PSPW WAVEFUNCTION_INITIALIZER

Listed below is the format of a WAVEFUNCTION_INITIALIZER sub-block.

PSPW
... 
   WAVEFUNCTION_INITIALIZER
     CELL_NAME: <string cell_name>
     WAVEFUNCTION_FILENAME: <string wavefunction_name default input_movecs>
     (RESTRICTED||UNRESTRICTED)
     if (RESTRICTED)   
        RESTRICTED_ELECTRONS: <integer restricted electrons>
     if (UNRESTRICTED) 
        UP_ELECTRONS: <integer up_electrons>
        DOWN_ELECTRONS: <integer down_electrons>
   END
...
END

The following list describes the input for the WAVEFUNCTION_INITIALIZER sub-block.

• $<$cell_name$>$ - name that points to the simulation_cell named $<$cell_name$>$. See section 31.2.3.
• $<$wavefunction_name$>$ - name that will point to a wavefunction file.
• RESTRICTED - keyword specifying that the calculation is restricted.
• UNRESTRICTED - keyword specifying that the calculation is unrestricted.

• $<$restricted_electrons$>$ - number of restricted electrons. Not used if an UNRESTRICTED calculation.
• $<$up_electrons$>$ - number of spin-up electrons. Not used if a RESTRICTED calculation.
• $<$down_electrons$>$ - number of spin-down electrons. Not used if a RESTRICTED calculation.

31.1.5.1 Old Style Input (version 3.3) to WAVEFUNCTION_INTITIALIZER

For backwards compatibility, the input to the WAVEFUNCTION_INITIALIZER sub-block can also be of the form

PSPW
... 
   WAVEFUNCTION_INITIALIZER
     CELL_NAME: <string cell_name>
     WAVEFUNCTION_FILENAME: <string wavefunction_name default input_movecs>
     (RESTRICTED||UNRESTRICTED)
     
     [UP_FILLING: <integer up_filling>
        [0 0 0   0]
        {<integer kx ky kz> (-2||-1||1||2)}]
     [DOWN_FILLING: <integer down_filling>
        [0 0 0   0]
        {<integer kx ky kz> (-2||-1||1||2)}]
   END
...
END

where

• $<$cell_name$>$ - name that points to the simulation_cell named $<$cell_name$>$. See section 31.2.3.
• $<$wavefunction_name$>$ - name that will point to a wavefunction file.
• RESTRICTED - keyword specifying that the calculation is restricted.
• UNRESTRICTED - keyword specifying that the calculation is unrestricted.
• $<$up_filling$>$ - number of restricted molecular orbitals if RESTRICTED and number of spin-up molecular orbitals if UNRESTRICTED.
• $<$down_filling$>$ - number of spin-down moleclar orbitals if UNRESTRICTED. Not used if a RESTRICTED calculation.
• $<$kx ky kz$>$ - specifies which planewave is to be filled.

The values for the planewave $(-2\vert\vert-1\vert\vert 1\vert\vert 2)$ are used to represent whether the specified planewave is a cosine or a sine function, in addition random noise can be added to these base functions. That is $+1$ represents a cosine function, and $-1$ represents a sine function. The $+2$ and $-2$ values are used to represent a cosine function with random components added and a sine function with random components added respectively.

31.1.6 V_WAVEFUNCTION_INITIALIZER

The v_wavefunction_initializer task is used to generate an initial velocity wavefunction datafile. Input to the V_WAVEFUNCTION_INITIALIZER task is contained within the V_WAVEFUNCTION_INITIALIZER sub-block.

PSPW
  ...
  V_WAVEFUNCTION_INITIALIZER
     ...
  END
  ...
END

To run a V_WAVEFUNCTION_INITIALIZER calculation the following directive is used:

TASK PSPW WAVEFUNCTION_INITIALIZER

Listed below is the format of a V_WAVEFUNCTION_INITIALIZER sub-block.

PSPW
... 
   V_WAVEFUNCTION_INITIALIZER
     V_WAVEFUNCTION_FILENAME: <string v_wavefunction_name default input_vmovecs>
     CELL_NAME: <string cell_name>
     (RESTRICTED||UNRESTRICTED)
     UP_FILLING: <integer up_filling>
     DOWN_FILLING: <integer down_filling>
   END
...
END

The following list describes the input for the V_WAVEFUNCTION_INITIALIZER sub-block.

• $<$cell_name$>$ - name that points to the simulation_cell named $<$cell_name$>$. See section 31.2.3.
• $<$wavefunction_name$>$ - name that will point to a velocity wavefunction file.
• RESTRICTED - keyword specifying that the calculation is restricted.
• UNRESTRICTED - keyword specifying that the calculation is unrestricted.
• $<$up_filling$>$ - number of restricted velocity molecular orbitals if RESTRICTED and number of spin-up velocity molecular orbitals if UNRESTRICTED.
• $<$down_filling$>$ - number of spin-down velocity moleclar orbitals if UNRESTRICTED. Not used if a RESTRICTED calculation.

31.1.7 WAVEFUNCTION_EXPANDER

The v_wavefunction_initializer task is used to convert a new wavefunction file that spans a larger grid space from an old wavefunction file. Input to the WAVEFUNCTION_EXPANDER task is contained within the WAVEFUNCTION_EXPANDER sub-block.

PSPW
  ...
  WAVEFUNCTION_EXPANDER
     ...
  END
  ...
END

To run a WAVEFUNCTION_EXPANDER calculation the following directive is used:

TASK PSPW WAVEFUNCTION_EXPANDER

Listed below is the format of a WAVEFUNCTION_EXPANDER sub-block.

PSPW
... 
   WAVEFUNCTION_EXPANDER   
     OLD_WAVEFUNCTION_FILENAME: <string old_wavefunction_name default input_movecs>
     NEW_WAVEFUNCTION_FILENAME: <string new_wavefunction_name default input_movecs>
     NEW_NGRID: <integer na1 na2 na3>
    
   END
...
END

The following list describes the input for the WAVEFUNCTION_EXPANDER sub-block.

• $<$old_wavefunction_name$>$ - name that points to a wavefunction file.
• $<$new_wavefunction_name$>$ - name that will point to a wavefunction file.
• $<$na1 na2 na3$>$ - number of grid points in each dimension for the new wavefunction file.

31.1.8 PSP_GENERATOR

A one-dimensional pseudopotential code has been integrated into NWChem. This code allows the user to modify and develop pseudopotentials. Currently, only the Hamann and Troullier-Martins norm-conserving pseudopotentials can be generated. This file can then be used by the pseudopotential_formatter task to generate a formatted pseudopotential file. Input to the PSP_GENERATOR task is contained within the PSP_GENERATOR sub-block.

PSPW
  ...
  PSP_GENERATOR
     ...
  END
  ...
END

To run a PSP_GENERATOR calculation the following directive is used:

TASK PSPW PSP_GENERATOR

Listed below is the format of a PSP_GENERATOR sub-block.

PSPW
... 
   PSP_GENERATOR
      PSEUDOPOTENTIAL_FILENAME: <string psp_name>
      ELEMENT: <string element>
      CHARGE: <real charge>
      MASS_NUMBER: <real mass_number>
      ATOMIC_FILLING: <integer ncore nvalence>
      ( (1||2||...) (s||p||d||f||...) <real filling> \
         ...)
      
      [CUTOFF: <integer lmax> 
         ( (s||p||d||f||g) <real rcut>\
         ...)
      ]
      PSEUDOPOTENTIAL_TYPE: (TROULLIER-MARTINS || HAMANN default HAMANN)
      SOLVER_TYPE: (PAULI || SCRHODINGER default PAULI)
      EXCHANGE_TYPE: (dirac || PBE96 default DIRAC)
      CORRELATION_TYPE: (VOSKO || PBE96 default VOSKO)
      [SEMICORE_RADIUS: <real rcore>]
      
   end
... 
END

The following list describes the input for the PSP_GENERATOR sub-block.

• $<$psp_name$>$ - name that points to a.
• $<$element$>$ - Atomic symbol.
• $<$charge$>$ - charge of the atom
• $<$mass$>$ - mass number for the atom
• $<$ncore$>$ - number of core states
• $<$nvalence$>$ - number of valence states.
• ATOMIC_FILLING:.....(see below)
• $<$filling$>$ - occupation of atomic state
• CUTOFF:....(see below)
• $<$rcore$>$ - value for the semicore radius (see below)

31.1.8.1 ATOMIC_FILLING Block

This required block is used to define the reference atom which is used to define the pseudopotential. After the ATOMIC_FILLING: $<$ncore$>$ $<$nvalence$>$ line The core states are listed (one per line), and then the valence states are listed (one per line). Each state contains two integer and a value. The first integer specifies the radial quantum number, $n$, The second integer specifies the angular momentum quantum number, $l$, and the third value specifies the occupation of the state.

For example to define a pseudopotential for the Neon atom in the $1s^2 2s^2 2p^6$ state could have the block

ATOMIC_FILLING: 1 2
        1  s  2.0   #core state    - 1s^2 
        2  s  2.0   #valence state - 2s^2
        2  p  6.0   #valence state - 2p^6

for a pseudopotential with a $2s$ and $2p$ valence electrons or the block

ATOMIC_FILLING: 3 0
        1  s  2.0    #core state
        2  s  2.0    #core state
        2  p  6.0    #core state

could be used for a pseudopotential with no valence electrons.

31.1.8.2 CUTOFF Block

This optional block specifies the cutoff distances used to match the all-electron atom to the pseudopotential atom. For Hamann pseudopotentials $r_{cut}(l)$ defines the distance where the all-electron potential is matched to the pseudopotential, and for Troullier-Martins pseudopotentials $r_{cut}(l)$ defines the distance where the all-electron orbital is matched to the pseudowavefunctions. Thus the definition of the radii depends on the type of pseudopotential. The cutoff radii used in Hamann pseudopotentials will be smaller than the cufoff raddi used in Troullier-Martins pseudopotentials.

For example to define a softened Hamann pseudopotential for Carbon would be

ATOMIC_FILLING: 1 2
  1  s  2.0
  2  s  2.0
  2  p  2.0
CUTOFF: 2
  s  0.8
  p  0.85
  d  0.85

while a similarly softened Troullier-Marting pseudopotential for Carbon would be

ATOMIC_FILLING: 1 2
  1  s  2.0
  2  s  2.0
  2  p  2.0
CUTOFF: 2
  s  1.200
  p  1.275
  d  1.275

31.1.8.3 SEMICORE_RADIUS Option

Specifying the SEMICORE_RADIUS option turns on the semicore correction approximation proposed by Louie et al (S.G. Louie, S. Froyen, and M.L. Cohen, Phys. Rev. B, 26, 1738, (1982)). This approximation is known to dramatically improve results for systems containing alkali and transition metal atoms.

The implmentation in the PSPW module defines the semi-core denisty, $\rho_{semicore}$ in terms of the core density, $\rho_{core}$, by using the sixth-order polynomial

$$\rho_{semicore}(r) = \left\{ \begin{array}{ll} \rho_{core} & \mbo... ...+ c_4 r^4 + c_5 r^5 + c_6 r^6 & \mbox{if $r < r_{semicore}$} \end{array}\right.$$ (31.1)

This expansion was suggested by Fuchs and Scheffler (M. Fuchs, and M. Scheffler, Comp. Phys. Comm.,119,67 (1999)), and is better behaved for taking derivatives (i.e. calculating ionic forces) than the expansion suggested by Louie et al.

31.2 PSPW RTDB Entries and DataFiles

Input to the PSPW module is contained in both the RTDB and datafiles. The RTDB is used to store input that the user will need to directly specify. Input of this kind includes ion positions, ion velocities, and simulation cell parameters. The datafiles are used to store input, such the one-electron orbitals, one-electron orbital velocities, formatted pseudopotentials, and one-dimensional pseudopotentials, that the user will in most cases run a program to generate.

31.2.1 Ion Positions

The positions of the ions are stored in the default geometry structure in the RTDB and must be specified using the GEOMETRY directive.

31.2.2 Ion Velocities

The velocities of the ions are stored in the default geometry structure in the RTDB, and must be specified using the GEOMETRY directive.

31.2.3 Simulation Cell

Simulation cells are stored in the RTDB. To enter a simulation cell into the RTDB the user defines a simulation_cell sub-block within the PSPW block. Listed below is the format of a simulation_cell sub-block.

PSPW
...
   SIMULATION_CELL
      CELL_NAME: <string name>
      BOUNDRY_CONDITIONS: (periodic || aperiodic)
      LATTICE_VECTORS:
        <real a1.x a1.y a1.z>
        <real a2.x a2.y a2.z>
        <real a3.x a3.y a3.z>
      NGRID: <integer na1 na2 na3>
   END
...
END

Basically, the user needs to enter the dimensions, gridding and boundry conditions of the simulation cell. The following list describes the input in detail.

• $<$name$>$ - user-supplied name for the simulation block.
• periodic - keyword specifying that the simulation cell has periodic boundary conditions.
• aperiodic - keyword specifying that the simulation cell has free-space boundary conditions.
• $<$a1.x a1.y a1.z$>$ - user-supplied values for the first lattice vector
• $<$a2.x a2.y a2.z$>$ - user-supplied values for the second lattice vector
• $<$a3.x a3.y a3.z$>$ - user-supplied values for the third lattice vector
• $<$na1 na2 na3$>$ - user-supplied values for discretization along lattice vector directions.

31.2.4 Analysis: Mulliken RTDB data

To perform Mulliken analysis information is needed from one-dimensional pseudopotential files. In-order to facilitate the transfer of this information to the simulation the ANALYSIS sub-block is used to exract the necessary information and put it into the RTDB.

PSPW
...
   ANALYSIS
      (psp_filename: <string psp_name> \
       ...)
   END
...
END

Basically, the user needs to enter each pseudopotential used in the simulation.

• $<$psp_name$>$ - name that ponts to a one-dimensional pseudopotential file.

31.2.5 Wavefunction Datafile

The one-electron orbitals are stored in a wavefunction datafile. This is a binary file and cannot be directly edited. This datafile is used by steepest_descent and Car-Parrinello tasks and can be generated using the wavefunction_initializer or wavefunction_expander tasks.

31.2.6 Velocity Wavefunction Datafile

The one-electron orbital velocities are stored in a velocity wavefunction datafile. This is a binary file and cannot be directly edited. This datafile is used by the Car-Parrinello task and can be generated using the v_wavefunction_initializer task.

31.2.7 Formatted Pseudopotential Datafile

The pseudopotentials in Kleinman-Bylander form expanded on a simulation cell (3d grid) are stored in a formatted pseudopotential datafile. This is a binary file and cannot be directly edited. This datafile is used by steepest_descent and Car-Parrinello tasks and can be generated using the pseudpotential_formatter task.

31.2.8 One-Dimensional Pseudopotential Datafile

The one-dimensional pseudopotentials are stored in a one-dimensional pseudopotential file. This is an ascii file and can be directly edited with a text editor. However, the user will usually use the psp_generator task to generate this datafile.

The data stored in the one-dimensional pseudopotential file is

   character*2 element       :: element name
   integer     charge        :: valence charge of ion
   real        mass          :: mass of ion
   integer     lmax          :: maximum angular component
   real        rcut(lmax)    :: cutoff radii used to define pseudopotentials
   integer     nr            :: number of points in the radial grid
   real        dr            :: linear spacing of the radial grid
   real        r(nr)         :: one-dimensional radial grid
   real        Vpsp(nr,lmax) :: one-dimensional pseudopotentials
   real        psi(nr,lmax)  :: one-dimensional pseudowavefunctions
   real        r_semicore        :: semicore radius
   real        rho_semicore(nr)  :: semicore density

and the format of it is:

[line 1:     ] element  
[line 2:     ] charge mass lmax
[line 3:     ] (rcut(l), l=1,lmax)
[line 4:     ] nr dr
[line 5:     ]    r(1)  (Vpsp(1,l),  l=1,lmax)
[line 6:     ] ....
[line nr+4:  ] r(nr) (Vpsp(nr,l), l=1,lmax)
[line nr+5:  ] r(1)  (psi(1,l), l=1,lmax) 
[line nr+6:  ] ....
[line 2*nr+4:] r(nr) (psi(nr,l), l=1,lmax)
[line 2*nr+5:] r_semicore
if (r_semicore read) then
[line 2*nr+6:] r(1)  rho_semicore(1)
[line 2*nr+7:] ....
[line 3*nr+5:] r(nr) rho_semicore(nr)
end if

31.2.9 PSPW Car-Parrinello Output Datafiles

31.2.9.1 XYZ motion file

Data file that stores ion positions and velocities as a function of time in XYZ format.

[line 1:          ]  n_ion
[line 2:          ]  
do ii=1,n_ion
[line 2+ii:       ] atom_name(ii), x(ii),y(ii),z(ii),vx(ii),vy(ii),vz(ii)
end do
[line n_ion+3     ] n_nion

do ii=1,n_ion
[line n_ion+3+ii: ] atom_name(ii), x(ii),y(ii),z(ii), vx(ii),vy(ii),vz(ii)
end do
[line 2*n_ion+4:  ]  ....

31.2.9.2 ION_MOTION motion file

Datafile that stores ion positions and velocities as a function of time

[line 1:          ]  it_out, n_ion, omega
[line 2:          ]  time
do ii=1,n_ion
[line 2+ii:       ] x(ii),y(ii),z(ii), vx(ii),vy(ii),vz(ii)
end do
[line n_ion+3     ] time
do 
do ii=1,n_ion
[line n_ion+3+ii: ] x(ii),y(ii),z(ii), vx(ii),vy(ii),vz(ii)
end do
[line 2*n_ion+4:  ]  ....

31.2.9.3 EMOTION motion file

Datafile that store energies as a function of time

[line 1:          ]  time, E1,E2,E3,E4,E5,E6,E7,E8, (E9,E10, if Nose-Hoover)
[line 2:          ] ...

31.2.9.4 HMOTION motion file

Datafile that stores the rotation matrix as a function of time.

[line 1:          ] time
[line 2:          ] ms,ne(ms),ne(ms)
do i=1,ne(ms)
[line 2+i:        ] (hml(i,j), j=1,ne(ms)
end do
[line 3+ne(ms):   ] time
[line 4+ne(ms):   ] ....

31.2.9.5 EIGMOTION motion file

Datafile that stores the eigenvalues for the one-electron orbitals as a function of time.

[line 1:          ]  time, (eig(i), i=1,number_orbitals)
[line 2:          ] ...

31.2.9.6 OMOTION motion file

Datafile that stores a reduced representation of the one-electron orbitals. To be used with a molecular orbital viewer that will be ported to NWChem in the near future.

31.3 Minimizing the DFT Energy Functional

In this section the algorithms for steepest_descent and conjugate_gradient methods that are used to minimize the DFT energy functional are presented. First is the the steepest descent method without line minimization. This method is slow, but doesn't require line minimization. Second is the recent method developed by Edelman, Arias, and Smith. These researchers have developed methods for performing line minimizations on the constrained space used in the DFT energy functional (A. Edelman, T. Arias, and S.T. Smith, Siam J. Matrix Anal. Appl., 20,303, (1998)).

31.3.1 Steepest Descent Equations

To minimize the DFT energy functional with respect to $\{\psi_{i,\sigma}\}$ we define the the auxiliary functional, $F$, which takes into account the orthonormality constraints.

$$F\left(\{\psi_{i,\sigma}\},\{\vec{R_I}\} \right) = E\left(\{\psi_{i,\sigma}\},\{\vec{R_I}\} \right)$$

$$+ \sum_{ij,\sigma} \left( \int d\vec{r}\ \psi_{i,\sigma}^{*}(\vec{r}) \psi_{j,\sigma}(\vec{r}) \Lambda_{ji,\sigma} - \delta_{ij,\sigma} \right)$$ (31.2)

The steepest descent minimization procedure without line minimization is then

$$\psi_{i,\sigma}^{t+ \Delta t} \leftarrow \psi_{i,\sigma}^{t} - \alpha \left[ \frac{\delta F}{\delta \psi_{... ...\psi_{i,\sigma}^{*}} + \sum_{j} \psi_{j,\sigma} \Lambda_{ji,\sigma} \right]_{t}$$ (31.3)

where $\alpha =\frac{\Delta t}{\sqrt{\mu}}$ is a positive numerical parameter. In this minimization procedure we have to know the variational derivative $\frac{\delta E}{\delta \psi_{i,\sigma}^{*}}$ and the matrix $\Lambda_{ij,\sigma}$. The variational derivative $\frac{\delta E}{\delta \psi_{i,\sigma}^{*}}$ can be analytically found and is

$$\frac{\delta E}{\delta \psi_{i,\sigma}^{*}} = -\frac{1}{2} \nabla^2 \psi_{i,\sigma}(\vec{r})$$

$$+ \int d\vec{r^{\prime}} W_{ext}(\vec{r},\vec{r^{\prime}}) \psi_{i,\sigma}(\vec{r^{\prime}})$$

$$+ \int d\vec{r^{\prime}} \frac{n(\vec{r^{\prime}})}{\vert\vec{r}-\vec{r^{\prime}}\vert} \psi_{i,\sigma}(\vec{r})$$

$$+ \mu_{xc}^{\sigma}(\vec{r}) \psi_{i,\sigma}(\vec{r})$$

$$\equiv \hat{H} \psi_{i,\sigma}$$ (31.4)

The ion positions can also be minimized using a fixed step minimization procedure.

$$\vec{R}_I^{t+\Delta t} \leftarrow \vec{R}_I^{t} - \frac{(\Delta t)}{\sqrt{M_I}} \frac{\partial E}{\partial \vec{R}_I}$$ (31.5)

To find the matrix $\Lambda_{ij,\sigma}$ we impose the orthonormality constraint on $\psi_{i,\sigma}^{t+\Delta t}$ to obtain the matrix Riccatti equation (to simplify the following equations we define the following symbol $\bar{\psi}_{i,\sigma}^{t} = \left(\frac{\delta E}{\delta \psi_{i,\sigma}^{*}} \right)_{t}$ )

$$I = <\psi_{i,\sigma}^{t+\Delta t} \vert \psi_{j,\sigma}^{t+\Delta t}>$$

$$= \left( <\psi_{i,\sigma}^{t}\vert - \alpha \left[ <\bar{\psi}_{i,\... ...rt + \sum_{k} \Lambda_{ik,\sigma}^{*} <\psi_{k,\sigma}^{t}\vert \right] \right)$$

$$\left( \vert\psi_{j,\sigma}^{t}> - \alpha \left[ \vert\bar{\psi}_... ...}^{t}> + \sum_{l} \vert\psi_{l,\sigma}^{t}> \Lambda_{lj,\sigma} \right] \right)$$

$$= <\psi_{i,\sigma}^{t}\vert \psi_{j,\sigma}^{t}>$$

$$- \alpha \left[ \begin{array}{ll} <\psi_{i,\sigma}^{t}\vert \bar{\p... ...{i,\sigma}^{t}\vert \psi_{l,\sigma}^{t}> \Lambda_{lj,\sigma} \end{array}\right]$$

$$+ \alpha^2 \left[ <\bar{\psi}_{i,\sigma}^{t}\vert \bar{\psi}_{j,\si... ...\bar{\psi}_{i,\sigma}^{t}\vert \psi_{l,\sigma}^{t}> \Lambda_{lj,\sigma} \right]$$

$$+ \alpha^2 \left[ \sum_{k} \sum_{l} \Lambda_{ik,\sigma}^{*} <\psi_{k,\sigma}^{t}\vert\psi_{l,\sigma}^{t}> \Lambda_{lj,\sigma} \right]$$

$$= A + \alpha XB + \alpha B^{\dag }X^{\dag } + \alpha^2 XCX^{\dag }$$ (31.6)

where $X_{ij,\sigma}=\Lambda_{ij,\sigma}^{*}$ and the matrices $A$, $B$, and $C$ are given by

$$A_{ij,\sigma} = \int d\vec{r} \left[\psi_{i,\sigma}^{t}(\vec{r}) - \alpha \left( ... ...) - \alpha \left( \frac{\delta E}{\delta \psi_{j,\sigma}^{*}} \right)_t \right]$$ (31.7)

$$B_{ij,\sigma} = \int d\vec{r} \left[\psi_{i,\sigma}^{t}(\vec{r}) \right]^{*} \lef... ...) - \alpha \left( \frac{\delta E}{\delta \psi_{j,\sigma}^{*}} \right)_t \right]$$ (31.8)

$$C_{ij,\sigma} = \int d\vec{r} \left[\psi_{i,\sigma}^{t}(\vec{r}) \right]^{*} \left[\psi_{j,\sigma}^{t}(\vec{r}) \right]$$ (31.9)

To solve the Riccatti equation an iterative solution is setup by rewriting Eq. 31.6 in the following form.

$$\alpha^2 XCX^{\dag } + \alpha X(B-I) + \alpha (B^{\dag }-I)X^{\dag } + \alpha(X + X^{\dag }) = I-A$$ (31.10)

Since X is by definition Hermitian, $X=X^{\dag }$, and

$$A \approx I + O(\alpha)$$

$$B \approx I + O(\alpha)$$ (31.11)

we can solve for X iteratively.

$$X_{(n)} \leftarrow \frac{1}{2 \alpha}(I-A) +$$

$$\frac{1}{2} \left( X_{(n-1)}(I-B) + (I-B^{\dag })X_{(n-1)}^{\dag } - \alpha X_{(n-1)}CX_{(n-1)}^{\dag } \right)$$ (31.12)

$$X_{(0)} \leftarrow \left(\frac{I-A}{2 \alpha}\right)$$ (31.13)

31.3.2 Conjugate Gradient with Curvature: Grassmann Manifold

Minimizing the DFT energy functional with a non-linear conjugate gradient algorithm requires an algorithm for minimizing along a search direction or line. However, the equation for a line (or search direction) in the DFT minimization problem is complicated by the fact that the energy functional is subject to orthonormality constraints. Furthermore, the orthonormality constraints complicate the conjugate gradient's algorithm for generating the current search direction by using a linear combination of the current gradient and the previous search direction, since the previous search direction must be parallel transported to the current search point.

These orthonormality constraint problems have been investigated and solved for by Edelman et al (A. Edelman, T. Arias, and S.T. Smith, Siam J. Matrix Anal. Appl., 20,303, (1998)). They recognized that the orthonormal Kohn-Sham orbitals $\{\psi_{i,\sigma}\}$ belong to a quotient space called a Grassmann manifold. Using the properties of the Grassmann manifold they developed an equation for a line (or geodesic) which conserves orthonormality, as well as a conjugate gradient procedure for generating the current search direction on this constrained space.

Following Edelman et al we write $\{ \psi_{i,\sigma}(\vec{r}) \}$ in terms of the orthonormal basis $\{ \phi_j(\vec{r}) \}$.

$$\psi_{i,\sigma}(\vec{r}) = \sum_{j}^{N_{basis}} \psi_{i,\sigma}(\phi_j) \phi_j(\vec{r})$$ (31.14)

Which allows us to rewrite the DFT energy functional as a function, $E(Y_{\uparrow},Y_{\downarrow})$, of two tall and skinny $N_{basis}$-by-$N_{\sigma}$ matrices which are

$$Y_{\uparrow}=\left[\begin{array}{rrrr} \psi_{1,\uparrow}(\p... ...{N_{\uparrow},\uparrow}(\phi_{N_{basis}}) \end{array} \right]$$ (31.15)

and

$$Y_{\downarrow}=\left[\begin{array}{rrrr} \psi_{1,\downarrow... ...\downarrow},\downarrow}(\phi_{N_{basis}}) \end{array} \right]$$ (31.16)

The orthonormality constraints make the $Y_{\sigma}$ matrices obey $Y_{\sigma}^{\dag }Y_{\sigma}=I$. Furthermore since the LSDA energy can be written as

$$E(Y_{\uparrow},Y_{\downarrow}) = \sum_{\sigma} tr \left( Y_{\sigma}^{\dag } (\mathrm{Operator}) Y_{\sigma} \right)$$ (31.17)

we are only interested in subspace spanned by $Y_\sigma$. This homogeneity property allows us to state that $E(Y_{\uparrow},Y_{\downarrow}) = E(Y_{\uparrow}Q_{\uparrow},Y_{\downarrow}Q_{\downarrow})$ where $Q_{\sigma}$ is any $N_{\sigma}$-by-$N_{\sigma}$ orthogonal matrix (i.e. $Q$ is any group element of the orthogonal group $O(N_{\sigma})$). The constrained surface that each $Y_{\sigma}$ spans is known as the Grassmann manifold (i.e. $Y_{\sigma}$ is a group element of the quotient space group $\frac{\left(\frac{O(N_{basis})}{O(N_{basis}-N_{\sigma})}\right)}{O(N_{\sigma})}$ ^31.1).

In this work Algorithm 31.1 (see below) is the conjugate gradient algorithm that is used to minimize the DFT energy functional $E=E\left( Y_{\uparrow},Y_{\downarrow} \right)$. What makes this algorithm different from a standard Polak-Ribière conjugate gradient algorithm is that a line search on a Euclidean space

$$Y_{\sigma}(t) = Y_{\sigma}^{(0)} + t*H_{\sigma}^{(0)}$$

is replaced by Eq. 31.19, and the parallel transports on a Euclidean space

$$\tau H_{\sigma}^{(0)} = H_{\sigma}^{(0)}$$

$$\tau G_{\sigma}^{(0)} = G_{\sigma}^{(0)}$$

are replaced by Eqs. 31.20-31.21.

Algorithm 31..1: Edelman et al's Algorithm for Conjugate Gradient Minimization on the Grassmann Manifold

Given $Y_{\sigma}^{(0)}$ such that $Y_{\sigma}^{(0)\dag }Y_{\sigma}^{(0)}$,
Compute $G_{\sigma}^{(0)} = \left[\frac{\delta E}{\delta Y_{\sigma}} \right]_{Y_{\sig... ...\left[\frac{\delta E}{\delta Y_{\sigma}} \right]_{Y_{\sigma}=Y_{\sigma}^{(0)}}$
Set $H_{\sigma}^{(0)} = -G_{\sigma}^{(0)}$. Find the compact singular value decompositions of $H_{\sigma}^{(0)} \rightarrow U_{\sigma} \Sigma_{\sigma} V_{\sigma}^{\dag }$ Minimize $E\left(Y_{\uparrow}(t), Y_{\downarrow}(t) \right)$ along the geodesic lines derived by Edelman et al for the Grassmann Manifold





Set and
compute Parallel transport along the geodesics the tangent vectors and







Compute the new search direction



where



Set , , and Go to step 2.

31.4 Car-Parrinello Scheme for Ab Initio Molecular Dynamics

Car and Parrinello developed a unified scheme for doing ab initio molecular dynamics by combining the motion of the ion cores and a fictacious motion for the Kohn-Sham orbitals of density-functional theory (R. Car and M. Parrinello, Phys. Rev. Lett. 55, 2471, (1985)). At the heart of this method they introduced a fictacious kinetic energy functional for the Kohn-Sham orbitals.

$$=$$ (31.21)

Given this kinetic energy the constrained equations of motion are found by taking the first variation of the auxiliary Lagrangian.

$$=$$ (31.22)

Which generates a dynamics for the wavefunctions and atoms positions through the constrained equations of motion:

$$=$$ (31.23)

$$=$$ (31.24)

where $\mu$ is the fictitious mass for the electronic degrees of freedom and are the ionic masses. The adjustable parameter $\mu$ is used to describe the relative rate at which the wavefunctions change with time. $\Lambda_{ij,\sigma}$ are the Lagrangian multipliers for the orthonormalization of the single-particle orbitals . They are defined by the orthonormalization constraint conditions and can be rigorously found. However, the equations of motion for the Lagrange multipliers depend on the specific algorithm used to integrate Eqs. 31.24-31.25.

For this method to give ionic motions that are physically meaningful the kinetic energy of the Kohn-Sham orbitals must be relatively small when compared to the kinetic energy of the ions. There are two ways where this criterion can fail. First, the numerical integrations for the Car-Parrinello equations of motion can often lead to large relative values of the kinetic energy of the Kohn-Sham orbitals relative to the kinetic energy of the ions. This kind of failure is easily fixed by requiring a more accurate numerical integration, i.e. use a smaller time step for the numerical integration. Second, during the motion of the system a the ions can be in locations where there is an Kohn-Sham orbital level crossing, i.e. the density-functional energy can have two states that are nearly degenerate. This kind of failure often occurs in the study of chemical reactions. This kind of failure is not easily fixed and requires the use of a more sophisticated density-functional energy that accounts for low-lying excited electronic states.

31.4.1 Verlet Algorithm for Integration

Eqs. 31.24-31.25 integrated using the Verlet algorithm results in

$$\psi_{i,\sigma}^{t+ \Delta t} \leftarrow$$ (31.25)

$$\vec{R}_I^{t+\Delta t} \leftarrow$$ (31.26)

In this molecular dynamic procedure we have to know variational derivative $\frac{\delta E}{\delta \psi_{i,\sigma}^{*}}$ and the matrix $\Lambda_{ij,\sigma}$. The variational derivative $\frac{\delta E}{\delta \psi_{i,\sigma}^{*}}$ can be analytically found and is

$$\frac{\delta E}{\delta \psi_{i,\sigma}^{*}} = -\frac{1}{2} \nabla^2 \psi_{i,\sigma}(\vec{r})$$

$$+ \int d\vec{r^{\prime}} W_{ext}(\vec{r},\vec{r^{\prime}}) \psi_{i,\sigma}(\vec{r^{\prime}})$$

$$+ \int d\vec{r^{\prime}} \frac{n(\vec{r^{\prime}})}{\vert\vec{r}-\vec{r^{\prime}}\vert} \psi_{i,\sigma}(\vec{r})$$

$$+ \mu_{xc}^{\sigma}(\vec{r}) \psi_{i,\sigma}(\vec{r})$$

$$\equiv \hat{H} \psi_{i,\sigma}$$ (31.27)

To find the matrix $\Lambda_{ij,\sigma}$ we impose the orthonormality constraint on $\psi_{i,\sigma}^{t+\Delta t}$ to obtain a matrix Riccatti equation, and then Riccatti equation is solved by an iterative solution (see section 31.3.1).

31.4.2 Constant Temperature Simulations: Nose-Hoover Thermostats

Nose-Hoover Thermostats for the electrons and ions can also be added to the Car-Parrinello simulation. In this type of simulation thermostats variables and are added to the simulation by adding the auxillary energy functionals to the total energy.

$$=$$ (31.28)

$$=$$ (31.29)

In these equations, the average kinetic energy for the ions is

(31.30)

where $f$ is the number of atomic degrees of freedom, is Boltzmans constant, and T is the desired temperature. Defining the average fictacious kinetic energy of the electrons is not as straighforward. Blöchl and Parrinello (P.E. Blöchl and M. Parrinello, Phys. Rev. B, 45, 9413, (1992)) have suggested the following formula for determining the average fictacious kinetic energy

(31.31)

where $\mu$ is the fictacious electronic mass, is average mass of one atom, and is the kinetic energy of the electrons.

Blöchl and Parrinello suggested that the choice of mass parameters, , and should be made such that the period of oscillating thermostats should be chosen larger than the typical time scale for the dynamical events of interest but shorter than the simulation time.

$$=$$ (31.32)

$$=$$ (31.33)

where $P_{ion}$ and $P_{electron}$ are the periods of oscillation for the ionic and ficatious electronic thermostats.

31.5 PSPW Tutorial 1: Minimizing the one-electron orbitals by Running a Steepest Descent and Conjuagate Gradient Simulation in Tandem

In this section we show how to setup and run a steepest descent simulation followed by a conjugate gradient simulation to optimize the one-electron orbitals with respect to the total energy for a NaCl molecule. It is practical to start most simulation in such a tandem fashion because the conjugate gradient optimizer does not work well when the one-electron orbitals are far from minimized.

Before running a steepest_descent or conjugate_gradient simulation several files and RTDB structures must be defined or initialized. Specifically, the user is required to have defined:

1. ion positions
2. simulation cell in the RTDB
3. formatted pseudopotential file for each kind of ion
4. file containing the one-electron orbitals
5. steepest_descent PSPW sub-block
6. conjugate_gradient PSPW sub-block

In the following tutorial we show the input needed to find an energy for an NaCl molecule. In this example default Hamann pseudopotentials are used for Na and Cl, the exchange correlation functional is LSDA, and the cutoff energy is 12 au. This example input deck can be found in the nwchem source tree in the file: nwchem/examples/pspw/NaCl.nw

1. Define the ion positions using the geometry directive.

            ...
            geometry units au 
            Na    2.23 0.0 0.0
            Cl   -2.23 0.0 0.0
            end
            ...
   

2. Generate a pseudopotential for Sodium atoms and store it in a one-dimensional pseudopotential datafile called Na.psp. The following input can be used to generate a default Hamann pseudopotential.

            ...
            PSPW
               PSP_GENERATOR
                  pseudopotential_filename: Na.psp
                   element: Na
                   charge: 11.0
                   mass_number: 23.0
                   solver_type: pauli
                   pseudopotential_type: hamann
                   atomic_filling: 3 1
                     1 s 2.0
                     2 s 2.0
                     2 p 6.0
                     3 s 1.0
               END
            END
            task PSPW PSP_GENERATOR
            ...
   

3. Use the same proceedure as above generate a pseudopotential for Chlorine.

            ...
            PSPW 
               PSP_GENERATOR
                  pseudopotential_filename: Cl.psp
                  element: Cl
                  charge: 17.0
                  mass_number: 35.0
                  solver_type: pauli
                  pseudopotential_type: hamann
                  atomic_filling: 3 2
                     1 s 2.0
                     2 s 2.0
                     2 p 6.0
                     3 s 2.0
                     3 p 5.0
               
               END
            END
            task PSPW PSP_GENERATOR
            ...
   

4. Define the simulation cell. The following input defines a simulation cell called ``small.'' This cell is periodic and cubic with a side length of 20.0 au and has 32 grid points in each direction.

            ...
            PSPW
               SIMULATION_CELL
                  cell_name: small
                  boundry_conditions: periodic
                  lattice_vectors:
                     20.0  0.0  0.0 
                      0.0 20.0  0.0 
                      0.0  0.0 20.0 
                  ngrid: 32 32 32
               END
            END
            ...
   

5. Format the Na pseudopotential onto the ``small'' simulation cell.

            ...
            PSPW
               PSP_FORMATTER
                  cell_name: small
                  psp_name: Na.psp
                  formatted_name: Na.vpp
               END
            END
            task PSPW PSP_FORMATTER
            ...
   

6. Format the Cl pseudopotential onto the ``small'' simulation cell.

            ...
            PSPW
               PSP_FORMATTER
                  cell_name: small
                  psp_name: Cl.psp
                  formatted_name: Cl.vpp
               END
            END
            task PSPW PSP_FORMATTER
            ...
   

7. Generate an initial guess for the one-electron orbitals.

            ...
            PSPW
               WAVEFUNCTION_INITIALIZER
                  cell_name: small
                  unrestricted
                  up_electrons: 4
                  down_electrons: 4
                  wavefunction_filename: NaCl.small.00.elc
               END
            END
            task PSPW WAVEFUNCTION_INTITALIZER
            ...
   

8. Do a coarse optimization of the one-electron orbitals with respect to energy using steepest descent.

            ...
            PSPW
               STEEPEST_DESCENT
                  cell_name: small
                  formatted_filename: Na.vpp
                  formatted_filename: Cl.vpp
                  input_wavefunction_filename:  NaCl.small.00.elc
                  output_wavefunction_filename: NaCl.small.00.elc
                  fake_mass: 400000.0d0
                  time_step: 51.8d0
                  loop: 10 100
                  tolerances: 1.0d-4 1.0d-4 1.0d-4
                  energy_cutoff:       21.0d0
                  wavefunction_cutoff: 21.0d0
               END 
            END
            task PSPW STEEPEST_DESCENT
            ...
   

9. Do a finer optimization of the one-electron orbitals with respect to energy using conjugate gradient, and perform a Mulliken analysis. Note that an analysis block must be defined.

            ...
            PSPW
              
               CONJUGATE_GRADIENT
                  cell_name: small
                  formatted_filename: Na.vpp
                  formatted_filename: Cl.vpp
                  input_wavefunction_filename:  NaCl.small.00.elc
                  output_wavefunction_filename: NaCl.small.00.elc
                  loop: 25 10
                  tolerances: 1.0d-9 1.0d-9 
                  energy_cutoff:       21.0d0
                  wavefunction_cutoff: 21.0d0
                  Mulliken
               END 
               ANALYSIS
                   psp_name: Na.psp
                   psp_name: Cl.psp
               END
            END
            task PSPW CONJUGATE_GRADIENT
            ...
   

31.6 PSPW Tutorial 2: Running a Car-Parrinello Simulation

In this section we show how to perform a Car-Parrinello molecular dynamic simulation for an NaCl molecule. As with the example in the previous tutorial, this example uses default Hamann pseudopotentials for Na and Cl, LSDA exchange correlation functional, and a cutoff energy of 12 au.

Before running a PSPW Car-Parrinello simulation the system should be on the Born-Oppenheimer surface, i.e. the one-electron orbitals should be minimized with respect to the total energy using steepest descent and conjugate gradient simulations(see section 31.5).

Continuning where we left off in section 31.5 we now continue the simulation to run a MD simulation of an NaCl molecule. This example input deck can be found in the nwchem source tree in the file: nwchem/examples/pspw/NaCl-md.nw

1. Optimize wavefunctions - See previous section

2. Generate an initial guess for the one-electron orbitals velocities.

            ...
            PSPW
               V_WAVEFUNCTION_INITIALIZER
                  cell_name: small
                  unrestricted
                  up_filling: 4
                  down_filling: 4
                  v_wavefunction_filename: NaCl.small.00.velc
               END
            END
            task PSPW WAVEFUNCTION_INTITALIZER
            ...
   

3. Run an MD simulation using Car-Parrinello.

            PSPW
               Car-Parrinello
                  cell_name: small
                  formatted_filename: Na.vpp
                  formatted_filename: Cl.vpp
                  input_wavefunction_filename:    NaCl.small.00.elc
                  output_wavefunction_filename:   NaCl.small.01.elc
                  v_input_wavefunction_filename:  NaCl.00.velc
                  v_output_wavefunction_filename: NaCl.01.velc
                  fake_mass: 800.0d0
                  time_step: 5.0d0
                  loop: 10 100
                  scaling: 1.0d0 1.0d0
                  energy_cutoff:       21.0d0
                  wavefunction_cutoff: 21.0d0
               END 
            END
            task PSPW Car-Parrinello
            ...
   

31.7 PSPW Capabilities and Limitations

• You cannot use more processors than the size of the third dimension (e.g. a 64x64x64 FFT grid can use at most 64 processors).
• The second and third dimensions of the FFT grid must be the same (i.e. the parameters na2 and na3 must be the same for each simulation cell).

31.8 Questions and Difficulties

Questions and encountered problems should be reported to nwchem-support@emsl.pnl.gov or to Eric J. Bylaska, Eric.Bylaska@pnl.gov