lp:~albertog/siesta/efield-2.5-jjunquera

Created by Alberto Garcia and last modified
Get this branch:
bzr branch lp:~albertog/siesta/efield-2.5-jjunquera
Only Alberto Garcia can upload to this branch. If you are Alberto Garcia please log in for upload directions.

Branch merges

Related bugs

Related blueprints

Branch information

Owner:
Alberto Garcia
Project:
Siesta
Status:
Development

Recent revisions

244. By Alberto Garcia

Implementation of forces with bulk electric fields

This commit updates the branch with a snapshot of the
state of Javier Junquera's code as of Feb 25, 2010.

A few cosmetic changes were needed for a simple compilation.
Note that the Examples directory has not been updated yet.

The detailed documentation for the implementation is not
yet committed.

added:
  Docs/efield.CHANGES
  Src/derdelk.f
  Src/grdelk.f
modified:
  Docs/CHANGES
  Src/Makefile
  Src/chempot.F
  Src/constr.f
  Src/delk.f
  Src/dfscf.f
  Src/dhscf.F
  Src/diagef.F
  Src/diagpol.f
  Src/efield.F
  Src/finitefield.F
  Src/indexkgrid.F
  Src/initdm.F
  Src/iodm.F
  Src/linpack.F
  Src/m_dimefield.f90
  Src/m_forces.F90
  Src/m_overkkisig.f
  Src/m_pkisigma.f
  Src/m_psitilde.f
  Src/meshmatrix.F
  Src/new_dm.F
  Src/overfsm.f
  Src/setup_hamiltonian.F
  Src/siesta.F
  Src/siesta_efield.F
  Src/siesta_forces.F
  Src/siesta_init.F
  Src/write_subs.F

243. By Alberto Garcia

Implementation of finite electric field for periodic systems (J Junquera)
--------------------------------
    In kgrid.F and find_kgrid.F
--------------------------------

     In order to have the k-points ordered,
 a vector of the reciprocal lattice has been subtracted to some of them.
 The algorithm is esentially the same, with the exception of the following
 lines:

C Find total range of grid indexes
C Find total range of grid indexes
! JJ: If the origin of the Monkhorst-Pack mesh has been displaced,
! substract a vector of the reciprocal lattice to some of the points,
! in order to have an ordered mesh of k-points.
! Nothing is done if the origin of the mesh of k-point is undisplaced.
      do j = 1,3
        ng(j) = kdsc(j,j)
        if (abs(dkg(j)) .lt. tiny) then
          igmin(j) = -( (ng(j)-1) / 2)
          igmax(j) = ng(j) / 2
        else
          igmin(j) = -( (ng(j)-1) / 2) - 1
          igmax(j) = ng(j) / 2 - 1
        endif
      enddo

      If the mesh of k-points is not displaced, nothing is done.
      If the mesh of k-points is displaced, and we have even numbers
 in the diagonal, the k-points are the same as before and only their order
 in the list is altered.
      If the mesh of k-points is displaced, and we have odd numbers
 in the diagonal, some k-points are different from the ones obtained
 with the original subroutine. Some of the components might have changed
 their sign.

--------------------------------
    In kpoint_grid.F90
--------------------------------
      The variables kscell and kdispl have been declared public
 (they are going to be passed to finiteef).

--------------------------------
    In alloc.F90:
--------------------------------

     New subroutines:
- realloc_z3
- realloc_z4
- dealloc_z3
- dealloc_z4

     Complex versions of the previous subroutines to reallocate or deallocate
 complex arrays of 3 and 4 dimensions.

     The complex pointers are supposed to be in every case
 of double precision type.

     realloc_z3 and realloc_z4 are copies of the previous
 realloc_d3 and realloc_d4 where the definition of the array has been changed
 from double precision to complex.

     dealloc_z3 and dealloc_z4 are copies of the previous
 dealloc_z1 where the dimension of the array to be deallocated
 has been increased to 3 and 4 respectively.
 from double precision to complex.

--------------------------------
    In dhscf.F:
--------------------------------

      New argument in the call to dhscf:

C integer iexpikr : Switch which fixes whether the matrix elements
C of exp(i k r) are calculated or not.
C 0 = not computed.
C 1 = compute the matrix elements of exp(i k r).
C -1 = compute the matrix elements of exp(-i k r).

      New call to the subroutine delk, where the matrix elements
 of a plane wave are computed in real space.

--------------------------------
    In meshsubs.F:
--------------------------------

    A new variable, iatfold, is defined in the module mesh.
 This is an integer variable of dimension 3 * number of atoms in the supercell.
 It is defined and computed independently of whether the finite field
 capability is used or not. The computational cost in CPU time and memory
 is negligeable, though.

    The variable is allocated and computed in the subroutine InitAtomMesh
 This variable stores the supercell vector required to fold the position
 of an atom back into the unit cell.

--------------------------------
    In siesta_options.F90
--------------------------------

    The default value for MD.TypeOfRun has been changed from Verlet to CG.

    Two new variables:

! logical :: bulkfinitefieldflag ! Determine whether there is a finite electric
! ! field in a periodic system or not.
! logical :: infinitefield ! Determines if the control of the program is
! ! taken by the module to compute the
! ! finite electric field or not.

--------------------------------
    In state_analysis.F
--------------------------------

     The KS energy is written only if the calculation is at zero
 external electric field in bulk.

     Under an external electric field, the electric enthalpy is
 written instead.

--------------------------------
    In local_DOS.F:
--------------------------------

     The calls to dhscf have been modified adding the
 new variable iexpikr in the list of arguments.

     Since the matrix elements of the plane wave are not required
 to compute the local DOS,
 iexpikr is set up to zero in local_DOS

--------------------------------
    In setup_hamiltonian.F
--------------------------------

     New variable:

     integer:: iexpikr ! Calculate matrix element of a plane wave exp(ikr)
                       ! 0 => no,
                       ! 1 => calculate exp(ikr)
                       ! -1 => calculate exp(-ikr)

     The calls to dhscf have been modified adding the
 new variable iexpikr in the list of arguments.

     Since the matrix elements of the plane wave are not required
 during the normal self-consistent-cicle,
 iexpikr is set up to zero in setup_hamiltonian.F

--------------------------------
    In siesta_analysis.F
--------------------------------

     The call to dhscf have been modified adding the
 new variable iexpikr in the list of arguments.

     Since the matrix elements of the plane wave are not required
 for the analysis,
 iexpikr is set up to zero in siesta_analysis.F

--------------------------------
    In grdsam.F:
--------------------------------

     The calls to dhscf have been modified adding the
 new variable iexpikr in the list of arguments.

     Since the matrix elements of the plane wave are not required
 during the normal self-consistent-cicle,
 iexpikr is always set up to zero in grdsam.F

--------------------------------
    In m_dimefield.f90:
--------------------------------

     Module where all the relevant variables,
 including:

     - Derived types for the list of k-point neighbours
 are defined.

--------------------------------
    In indexkgrid.F
--------------------------------

      New subroutine that initializes some lists that allow us
 to relate the index of a given k-point with the indices
 of neighbouring k-points along the three directions of the
 reciprocal lattice.

      Look for instance at Eq. (53) of the paper by
 I. Souza et al., Phys. Rev. B 69, 085106 (2004),
 where the overlap matrices between a given k-point and
 the six first neighbours are required.

      The algorithm is inspired in the Appendix
 of the technical paper of Siesta.
 Basically:

- We reintroduce again the points that were cutted due to the
  time reversal symmetry (k = -k).
- Construct a supercell in reciprocal space.
- Look for neighboring points in the supercell.
- Once the neighbouring point is localized, go back to the normal mesh
 of k-points.
- Localize the equivalent point due to time reversal symmetry.

      Parallelization issues:
 The subroutine is called only once, and it takes no more than one second
 to compute the list of k-point neighbours.
 Therefore, no parallelization is required.

 The only place where parallelization might play a role
 is in the broadcasting of the derived types
 dictionaryneighkafter and dictionaryneighkbefore to all the nodes.

--------------------------------
    In delk.f
--------------------------------

      New subroutine where the matrix elements of a plane wave are
 computed in real space.

      The matrix elements of a plane wave in a basis of numerical
 atomic orbitals are computed in the grid, instead of in
 the former approach implemented by Daniel Sanchez-Portal,
 where a Taylor expansion of the plane wave was made assuming
 small values of the exponent:

      exp(i k r) \sim 1 + i k r

      The reason for this change is the fact that, within the new
 algorithm, there is no reason to assume that k (or more precisely
 the vector between consecutive k-points \Delta k)
 would be small enough to justify the cut of the Taylor expansion.
 Remember that, in the algorithm described by I. Souza et al.,
 the larger the number of k-points (that is, the smaller the
 interval between neighboring points) the smaller the value of the
 threshold energy for Zener tunneling
 [see the end of the first paragraph in the
 left column of the third page of the paper by
 I. Souza et al., Phys. Rev. Lett. 89, 117602 (2002)].

      Parallelization issues:
 If vmat is parallelized, delk should work as well, since the
 basic loops are maintained without any change.

--------------------------------
    In efield.F
--------------------------------

      New variables:

! logical bulkfinitefieldflag: switch that determines whether we are going to
! perform a calculation of the bulk system
! under the presence of an external electric field.

      Behaviour:

      Read the electric field and generate the corresponding potential
 in the grid if and only if we are not dealing with periodic systems.

      If we are dealing with an open system
 (molecule, cluster, wire, tube, surface)
 and we have a field normal to the periodic direction, then
 the potential corresponding to the field is introduced as a ramp with
 the discontinuity in the vacuum, in the usual way.

--------------------------------
    In m_energies.F90
--------------------------------

     New variable: DUextbulk (interaction energy with external electric field
                   in a bulk material)

--------------------------------
    In post_scf_work.F
--------------------------------

     New term added to the sum of the total energy:

      Etot = E0 + DEna + DUscf + DUext + Exc + Ecorrec + Emad + Emm +
     . Emeta + DUextbulk

--------------------------------
    In scfconvergence_test.F
--------------------------------

     New term added to the sum of the total energy:

      Etot = E0 + DEna + DUscf + DUext + Exc + Ecorrec + Emad + Emm +
     . Emeta + DUextbulk

--------------------------------
    In siesta_init.F
--------------------------------

     The variable DUextbulk is initialized

     DUextbulk = 0.0_dp

--------------------------------
    In write_subs.F
--------------------------------

     The variable DUextbulk is written in the same way as DUext was
 written before.

     Besides, the writing of the decomposition of the energy,
 forces, and stresses has been changed. Now, the "prefix" of everyline
 (as it was now "siesta :") depends on whether the calculation
 is at zero electric field (in this case "siesta :") is maintained,
 or at a finite electric field (in this case, the prefix is "finief:").

     The Harris energy is only written under zero electric field in bulk.

--------------------------------
    m_noccbands.f
--------------------------------

      New module (extracted from Daniel's subroutine ksv.f)
 that computes the number of occupied bands.

--------------------------------
    diagef.F
--------------------------------

      New subroutine of diagonalization.

      All the matrices are defined as complex arrays.
      The matrix that is diagonalized is Tk and not the KS hamiltonian.
      It builds the density and energy density matrices,
 but now the loop in k-points is outside the subroutine
 inside finitefield).

--------------------------------
    m_overkisig.f
--------------------------------

      File that contains the subroutine that computes the overlap
 between the periodic part of the Bloch-like resonances in neighboring
 k-points, following the definition given right before
 Eq. (53) of the paper by I. Souza et al. Phys. Rev. B 69, 085106 (2004).

      To adapt this to Siesta, we have followed Eq. (5) of the paper
 by D. Sanchez-Portal et al., Fundamental Physics of Ferroelectrics,
 535, 111, 2000.

--------------------------------
    m_psitilde.f
--------------------------------

     Compute the dual basis of the Bloch-like resonances,
 vtilde [see Eq. (53) of the paper by
 Ivo Souza et al. Phys. Rev. B 69, 085106, (2004) ]

--------------------------------
    m_pkisigma.f
--------------------------------

     Compute the projectors defined in Eq. (61) of the paper by
 Ivo Souza et al. Phys. Rev. B 69, 085106, (2004)

--------------------------------
    siesta_efield.F
--------------------------------

      Interface between siesta.F and the module that computes
 the response of a dielectric to a finite electric field.

242. By Alberto Garcia

tag of <email address hidden>/siesta-devel--reference--2.5--patch-11

241. By Alberto Garcia

Proper fix of cell management for zmatrix case. Enthalpy cosmetics
A recent previous patch broke the unit cell handling in the zmatrix case.
It has been fixed.

The assorted "enthalpy" output has been rationalized,
reserving the use of the word "enthalpy" for the original case
of "target enthalpy".

(Fixed syntax error in std-test.mk)

240. By Alberto Garcia

Make diagonalization the default solution method regardless of size
Make diagonalization the default solution method regardless of system size.

239. By Alberto Garcia

Fix test.mk scripts
The std-test.mk script did not copy the output file to the main test directory.
The bsc-test.mk did not remove all the scratch files.

Fixed a typo in the input to the si2x1h test.

238. By Alberto Garcia

Fix for structural supercell construction
Due to the accidental mis-use of a module variable, the routine coor
did not generate the supercell correctly when a %supercell block was used.

Added and additional warning if such a block is used when the structure is
input in Zmatrix form.

237. By Alberto Garcia

Fix memory leak in dhscf
The array dvxcdn was allocated but never deallocated in dhscf. Actually,
it was never used in cellxc, so it has now been dimensioned trivially
to (1,1,1). The bug was introduced in the bsc-master-2.1 branch, at
patch 27.

236. By Alberto Garcia

New 'fdf' target in Makefile for tests.
(E. Anglada)

New target for the test option: fdf calls.
Requested by BSC so they can test the new fdf.

New files (new tests and missing xml output files for reference)
A/ Tests/var_cell_stress
A Tests/Reference-xml/fire_benzene.xml
A Tests/Reference-xml/si_coop.xml
A Tests/var_cell_stress/.arch-ids/=id
A Tests/var_cell_stress/.arch-ids/makefile.id
A Tests/var_cell_stress/.arch-ids/var_cell_stress.fdf.id
A Tests/var_cell_stress/.arch-ids/var_cell_stress.pseudos.id
A Tests/var_cell_stress/makefile
A Tests/var_cell_stress/var_cell_stress.fdf
A Tests/var_cell_stress/var_cell_stress.pseudos

Modified files
M Src/siesta_move.F : remove debug info
M Src/remove_intramol_pressure.f90 : remove debug info
M Tests/std-Makefile Xmlcheck enabled by default. New targets: fdf, all
M Tests/std-test.mk
M Tests/si2x1h/si2x1h.fdf New option
M Src/atom.f Remove debug info.

(replayed patch 6 from outlier branch ref-2.4)

235. By Alberto Garcia

xml tester update (Eduardo Anglada)
New files:
A Src/xmlparser/corresponding_node.f90
A Src/xmlparser/string_utilities.f90

Changes in the xml output:
M Src/state_analysis.F
M Src/siesta_move.F
M Src/setup_hamiltonian.F
M Src/siesta_forces.F
M Src/scfconvergence_test.F
M Src/post_scf_work.F

Changes in the xml tester:
M Src/wxml/m_wxml_core.f90
M Src/wxml/flib_wstml.f90
M Src/xmlparser/makefile
M Src/xmlparser/m_dom_element.f90
M Src/xmlparser/m_strings.f90
M Src/xmlparser/compare_m.f90
M Src/xmlparser/test.f90

Changes in the tests makefiles
M Tests/std-Makefile
M Tests/std-test.mk

Upgraded xml reference outputs:
...

(replayed patch-5 from outlier branch ref-2.4)

Branch metadata

Branch format:
Branch format 7
Repository format:
Bazaar repository format 2a (needs bzr 1.16 or later)
Stacked on:
lp:siesta
This branch contains Public information 
Everyone can see this information.

Subscribers