lp:~nickpapior/siesta/low-D-stress

Branch merges

Propose for merging

No branches dependent on this one.

Ready for review for merging into lp:siesta/4.1

Alberto Garcia: Pending requested 2019-02-14

Diff: 2798 lines (+2016/-685) (has conflicts)

12 files modified

Src/Makefile (+5/-3)
Src/write_subs.F90 (+884/-681)
Tests/Makefile (+1/-1)
Tests/Reference/graphene_non_perp.out (+521/-0)
Tests/Reference/graphene_perp.out (+521/-0)
Tests/graphene_non_perp/graphene_non_perp.fdf (+33/-0)
Tests/graphene_non_perp/graphene_non_perp.pseudos (+1/-0)
Tests/graphene_non_perp/makefile (+6/-0)
Tests/graphene_perp/graphene_perp.fdf (+33/-0)
Tests/graphene_perp/graphene_perp.pseudos (+1/-0)
Tests/graphene_perp/makefile (+6/-0)
version.info (+4/-0)

api

Related source package recipes

No recipes using this branch. (?)

Create packaging recipe

Related snap packages

Related bugs

Related blueprints

1072. By Nick Papior on 2019-02-28

Finalized calculation of low-D stress tensors using the shaper routine

We now calculate the vacuum regions based on the shaper routine.
This results in very similar output for skewed lattice vectors.
See e.g. the two new tests.

  graphene_perp
  graphene_non_perp

1071. By Nick Papior on 2019-02-26

Added two examples where this should be printed out

Currently I have not added the output since it may change.

1070. By Nick Papior on 2019-02-14

Added 1D and 2D stress tensor calculations

The problem with low-D materials is that stress is "per volume",
but since vacuum is independent in a slab calculation the stress will
be variable depending on the amount of vacuum.

What one *really* should do is to use the volume of the slab along the
vacuum direction so it may be fixed depending on the extend of the slab.
This poses additional problems since how do you decide the slab thickness?
- vdW bond-lengths
- orbital radii
- or 3rd?

This amends output to do the simplest thing but forces users to do post-processing.
The vacuum region is divided out and thus change the units depending on the dimensionality.
- 2D == eV/Ang^2
- 1D == eV/Ang

Then the user can them-selves determine the actual height/area of the slab/chain.

Note that the regular stress tensors are still printed out, but now an additional
block is added.

TODO: these stress's should probably be the ones used for relaxations and var-cell
      calculations?

1069. By Nick Papior on 2019-02-14

Cleaned up write_subs.F90 and changed implicit loops to explicit notations

This merely makes the code more readable and makes it ready for the low-D changes.

Fixed all use statements in write_subs.
Moved siesta_write_forces to a separate module, along the lines of the pressure,
energetics, etc.

1068. By Nick Papior on 2019-02-14

Moved write_subs.F to write_subs.F90 for easier maintenance

1067. By Nick Papior on 2019-02-13

Added a new Bloch-expansion algorithm

It does exactly the same thing but is somewhere between 1 and 2 orders of
magnitudes faster.

The reason is that the Bloch-expansion has a Toeplitz structure allowing for
reduced calculations and copying data around. This is now being done at
two levels:

1) For Bloch expansions along 1 dimension we have an efficient algorithm.
2) For Bloch expansions along 2 dimensions we have an efficient algorithm.

Both are also better parallelized using OpenMP. Even for simple Bloch expanions
of ~5 I see a halving of execution time of the methods. This could potentially be
very important for large scale TB calculations, but also for regular TS calculations.

This also adds warnings when using the repetition of bloch expansion since
the algorithm only works for tiling constructs. Hence, when going from repeating
to tiling one would find an even larger performance improvement.

Added a better error message when the electrode coordinates are not coinciding.
A direct comparison between the device and electrode coordinates are
now being done.

1066. By Alberto Garcia on 2019-02-12

Simplify the choice between SiestaXC and BSC's versions of cellxc

The BSC_CELLXC preprocessor blocks have been removed and replaced
by conditional blocks which can be selected at runtime.

The user can control which version of cellxc is used by means of the
fdf logical variable:

     XC.Use.BSC.CellXC (T/F)

which is 'false' by default, thus using SiestaXC's cellxc.

The BSC version is slightly more efficient for GGA operations, and the
SiestaXC version must (and will) be used when dealing with van der
Waals density functionals.

NOTES:

* BSC's version of cellxc now uses ldaxc, ggaxc, and setxc/getxc from
  SiestaXC. This enhances its functionality (more functionals) and
  saves on duplicated code. The old 'xc.f' file has been removed.

* In 'meshsubs', 'distriphionmesh' now accepts an extra argument to
  flag the need for stencil initilization.

* SiestaXC is now seeded with the "linear" pre-defined parallel
  distribution instead of having it use its internal load balancer.

* In 'forhar', a number of 'reord' operations towards 'sequential'
  fine-point ordering (which were compiled-in only for the new
  interface) have been removed, as all the arrays involved are already
  in 'sequential' form.

  Tests show that these changes do not seem to affect the results, however.

  The 'forhar' interface uses the 'ntm' argument in all cases, to simplify
  the code.

* Removed bogus incompatibility check in ldau.F

1065. By Alberto Garcia on 2019-02-12

Merge 4.0:591, with added flexibility in the handling of pseudopotentials

* The maximum l for KB projectors is now automatically set to the
highest l in the pseudopotential file. The program will use Vlocal to
construct the PAOs if there are not enough l-channels of appropriate
semilocal pseudopotentials.

* The reading of the PS.KB-projectors block has been made more robust,
actually allowing missing reference-energy lines. 'nkbmx' has been
increased to 3 in 'atmparams'.

(+ Fix for remaining problem in assignment of n-quantum number to
polarization orbitals. Removal of dead code.)

(+ Fix for duplicate symbols in Util/Gen-basis)

1064. By Nick Papior on 2019-02-06

Clarified and fixed the complex parts of non-equilibrium contours

When using different eta values in electrodes the TSCCNEQ-* files contained
wrong imaginary parts. This commit changes the TSCCNEQ files to the following
naming scheme:

  <syslabel>.TSCCNEQ-<ID>-<mu>_<electrode>

meaning that the ID is the *internal* ID number and is useful for determining
the order of the energy-points.
The <mu> denotes that the contour is the non-equilibrium corrections to the
chemical potential named <mu>.
The <electrode> denotes that the weights in this file are associated
with electrode <electrode> and the imaginary part corresponds to the
imaginary part used in the self-energy calculation.
Previously all of these files contained the same imaginary part equal to the
TS.Contours.nEq.Eta value (and thus had no relation to the actually used eta
value for the electrode).

When the eta values are the same in all electrodes, the files will
contain the same values, however, when eta differs, so does the *imaginary*
part in the corresponding *_<electrode> files.

Conversion factors has been revised in the contour parts since it was
difficult to get the equivalent precision using external tools.
This is because the fdf-flags read in (which determines the energies)
has one conversion factor (eV->Ry) while the written files was using
the conversion factor in units.f90. This meant that there was a relative
difference of 1e-6 between the two. Now the contours are written using
the fdf-conversion factor, explicitly.

The precision in the written contour files has been increased to
maximum numerical precision. width == 17.

Added a check that all non-equilibrium electrode Eta values are not 0.

NOTE: I have asked Alberto what to do about unit conventions in the code.

1063. By Nick Papior on 2019-02-05

Updated install_* library versions, netcdf and hdf5.