Merge branch 'tk-thread-calc-matrix' of github.com:OrderN/CONQUEST-release into tk-thread-calc-matrix
This commit is contained in:
Tuomas Koskela
commit
d75229771c
|
|
@ -0,0 +1,518 @@
|
|||
.. _ase-conquest:
|
||||
|
||||
==========================================
|
||||
Managing Conquest with ASE
|
||||
==========================================
|
||||
|
||||
Below we give an introduction how to setup the ASE environment with respect
|
||||
to CONQUEST repository along with a few examples of ASE/Conquest capabilities.
|
||||
We assume that a python script or jupyter-notebook is used.
|
||||
|
||||
.. _ase_setup:
|
||||
|
||||
Setup
|
||||
-----
|
||||
|
||||
.. _ase_env:
|
||||
|
||||
Environment variables
|
||||
+++++++++++++++++++++
|
||||
|
||||
The script will need to set environmental variables specifying the
|
||||
locations of the CONQUEST executable ``Conquest``, and if required, the basis
|
||||
set generation executable ``MakeIonFiles`` and pseudopotential database.
|
||||
These variables are:
|
||||
|
||||
* ``ASE_CONQUEST_COMMAND``: the Conquest executable command including MPI/openMPI prefix.
|
||||
* ``CQ_PP_PATH``: the PAO path directory to where are located the the ``.ion`` files.
|
||||
* (optional) ``CQ_GEN_BASIS_CMD`` : the PAO generation executable ``MakeIonFiles``.
|
||||
|
||||
..
|
||||
Given the Conquest root directory ``CQ_ROOT`` which contains
|
||||
|
||||
..
|
||||
::
|
||||
|
||||
bin/
|
||||
docs/
|
||||
pseudo-and-pao/
|
||||
src/
|
||||
testsuite/
|
||||
tools/
|
||||
|
||||
Given the Conquest root directory ``CQ_ROOT``, initialisation might look to something like
|
||||
::
|
||||
|
||||
import os
|
||||
|
||||
CQ_ROOT = 'PATH_TO_CONQUEST_DIRECTORY'
|
||||
|
||||
os.environ['ASE_CONQUEST_COMMAND'] = 'mpirun -np 4 '+CQ_ROOT+'/bin/Conquest'
|
||||
os.environ["CQ_GEN_BASIS_CMD"] = CQ_ROOT+'/bin/MakeIonFiles"
|
||||
os.environ['CQ_PP_PATH'] = CQ_ROOT+'/pseudo-and-pao/'
|
||||
|
||||
|
||||
Go to :ref:`top <ase-conquest>`.
|
||||
|
||||
.. _ase_pao:
|
||||
|
||||
Pseudopotential/PAO files
|
||||
+++++++++++++++++++++++++
|
||||
|
||||
Conquest atomic pseudotential and basis functions are store in the ``.ion``
|
||||
files which will ne referred as to PAO files. Provided the pseudopotential files ``.pot`` available in ``CQ_PP_PATH``,
|
||||
automatic generation of numerical PAOs is possible using the program :ref:`MakeIonFiles <generating_paos>`
|
||||
available from the Conquest package.
|
||||
|
||||
Provided the PAO files, the basis set is specified through a python dictionary,
|
||||
for example:
|
||||
::
|
||||
|
||||
basis = {'O' : {'file': 'O_SZP.ion'},
|
||||
'H' : {'file': 'H_SZP.ion'},
|
||||
'C' : {'file': 'C_SZP.ion'}}
|
||||
|
||||
In this case they are all assumed to be obtained from Hamann pseudopotentials,
|
||||
which are the default. Knowing the the exchange and correlation functional <XC>
|
||||
from the Conquest input (*vide infra*) and the chemical symbol <X>, the *Calculator*
|
||||
will search the ``.ion`` file in different places:
|
||||
::
|
||||
|
||||
CQ_PP_PATH
|
||||
CQ_PP_PATH/lib/
|
||||
CQ_PP_PATH/<XC>/<X>/
|
||||
|
||||
including the current directory and the ASE working directory (*vide infra*). If
|
||||
your PAO file is located in a different place you can include the path in the
|
||||
basis dictionary:
|
||||
::
|
||||
|
||||
basis = {'O' : {'file': 'O_SZP.ion',
|
||||
'directory': '<PATH_TO_FILE>'},
|
||||
'H' : {'file' : 'H_SZP.ion'},
|
||||
'C' : {'file' : 'C_SZP.ion'}}
|
||||
|
||||
|
||||
For generating the PAO files, the keyword ``gen_basis`` should be set to ``True``
|
||||
(default is ``False``) and the size be provided (default is ``medium``).
|
||||
For instance:
|
||||
::
|
||||
|
||||
basis = {'O' : {'gen_basis' : True,
|
||||
'basis_size': 'small'},
|
||||
'C' : {'gen_basis' : True,
|
||||
'basis_size': 'medium'},
|
||||
'H' : {'gen_basis' : True,
|
||||
'basis_size': 'large'}}
|
||||
|
||||
will create the ``O.ion``, ``C.ion`` and ``H.ion`` files where ``small``, ``medium`` and ``large``
|
||||
are :ref:`default size basis set <pao_gen_default>`. You are allowed to choose the functional and
|
||||
to add options for basis set generation:
|
||||
::
|
||||
|
||||
basis = {'H' : {'gen_basis' : True,
|
||||
'basis_size': 'small',
|
||||
'xc' : 'LDA',
|
||||
'Atom.Perturbative_Polarised': False}}
|
||||
|
||||
.. note::
|
||||
Only Hamann pseudopotentials for LDA, PBE and PBEsol are available within
|
||||
the CONQUEST distribution. For using other functionals see
|
||||
:ref:`Generating new pseudopotentials <pao_gen_pseudo>`.
|
||||
|
||||
.. warning::
|
||||
Generating polarised PAOs for some atoms can be problematic (mainly group I
|
||||
and II). Please review carefuly the ``MakeIonFiles`` input files
|
||||
named ``Conquest_ion_input`` which are collected in ``CQ_PP_PATH/<XC>/<X>/``
|
||||
if you are not sure about what you are doing, and check your PAOs.
|
||||
|
||||
|
||||
Go to :ref:`top <ase-conquest>`.
|
||||
|
||||
.. _ase_calculator:
|
||||
|
||||
CONQUEST Calculator
|
||||
--------------------
|
||||
|
||||
The CONQUEST *Calculator* class can be invoked from the ase Calculator set as described
|
||||
in the example below:
|
||||
::
|
||||
|
||||
from ase.calculators.conquest import Conquest
|
||||
|
||||
A minimal example is given below for setting the CONQUEST *Calculator* (named ``calc``)
|
||||
of the `ASE Atoms object <https://wiki.fysik.dtu.dk/ase/ase/atoms.html#module-ase.atoms>`_
|
||||
named `struct`:
|
||||
::
|
||||
|
||||
from ase.calculators.conquest import Conquest
|
||||
from ase.build import bulk
|
||||
|
||||
struct = bulk('NaCl', crystalstructure='rocksalt', a=5.71, cubic=True)
|
||||
basis = {'Cl' : {'file' : 'Cl.ion'}, 'Na' : {'file' : 'Na.ion'}}
|
||||
|
||||
calc = Conquest(basis=basis,atoms=struct)
|
||||
|
||||
or, equivalently,
|
||||
::
|
||||
|
||||
from ase.calculators.conquest import Conquest
|
||||
from ase.build import bulk
|
||||
|
||||
struct = bulk('NaCl', crystalstructure='rocksalt', a=5.71, cubic=True)
|
||||
basis = {'Cl' : {'file' : 'Cl.ion'}, 'Na' : {'file' : 'Na.ion'}}
|
||||
|
||||
struct.calc = Conquest(basis=basis)
|
||||
|
||||
|
||||
In basic calculate mode (compute ``energy``), the *Calculator* comes with 3 *methods*:
|
||||
|
||||
* ``write_input()``:
|
||||
this function will setup the :ref:`input files <setting_up>`. For CONQUEST, the PAO basis will be generated/copied with respect to the dictionary key/value pairs, and ``Conquest_input`` file including the calculation parameters will be written, a long with the ``coordinate`` file, containing the lattice vectors (in Bohr Unit) and atomic positions (in fractional coordinates).
|
||||
|
||||
* ``execute()``:
|
||||
this function execute the calculation. For CONQUEST, it will launch the ``ASE_CONQUEST_COMMAND`` setup in :ref:`the environment varaibles <ase_env>`.
|
||||
|
||||
* ``read_results()``:
|
||||
this function post-process the the output file. For CONQUEST, the ``energy``, ``forces``, ``stress`` and ``eigenvalues`` will be extracted from the ``Conquest_out_ase`` output file.
|
||||
|
||||
.. note::
|
||||
The funtion ``read_results()`` operate on the ``Conquest_out_ase`` file.
|
||||
This output file is not created by default by CONQUEST. If you want to post-process
|
||||
a calculation with an input generated by hand you must add ``IO.WriteOutToASEFile True``
|
||||
in ``conquest_input``.
|
||||
|
||||
|
||||
The **indirect** way for managing CONQUEST calculation with ASE is:
|
||||
::
|
||||
|
||||
struct.calc = Conquest(basis=basis)
|
||||
|
||||
struct.calc.write_input(struct)
|
||||
struct.calc.execute()
|
||||
struct.calc.read_results(struct)
|
||||
|
||||
|
||||
where ``struct.calc.execute()`` can be ignored when, for instance, the calculation
|
||||
is performed on a supercomputer and the output file is then copied back to the current
|
||||
directory for post-processing.
|
||||
|
||||
|
||||
The **direct** way is simply:
|
||||
::
|
||||
|
||||
struct.calc = Conquest(basis=basis)
|
||||
|
||||
struct.calc.calculate(struct)
|
||||
|
||||
or, equivalently,
|
||||
::
|
||||
|
||||
struct.calc = Conquest(basis=basis)
|
||||
|
||||
struct.get_potential_energy()
|
||||
|
||||
Go to :ref:`top <ase-conquest>`.
|
||||
|
||||
.. _ase_input:
|
||||
|
||||
Keywords for generating the Conquest_input file
|
||||
-----------------------------------------------
|
||||
|
||||
In principle all the `Conquest input parameters <https://conquest.readthedocs.io/en/latest/input_tags.html>`_
|
||||
can be added to ``Conquest_out_ase`` using key/value pairs in a dictionary. There are 3 class of parameters:
|
||||
|
||||
* mandatory : they are parsed to the *Calculator* and have no defaults ; there are mandatory.
|
||||
* important : they are parsed to the *Calculator* they can be freely modified. Some of them are pure ASE keywords.
|
||||
* defaults : they are set as defaults ; some of them must not be modified. They are read by the *Calculator* through a dictionay ``conquest_flags``.
|
||||
|
||||
|
||||
Mandatory keywords
|
||||
++++++++++++++++++
|
||||
|
||||
|
||||
=================== ===================== =============== ================================
|
||||
keyword type default value description
|
||||
=================== ===================== =============== ================================
|
||||
``atoms`` ``atoms`` None an atoms object constructed either via ASE or read from an input
|
||||
``basis`` ``dict`` None a dictionary specifying the pseudopotential/basis files
|
||||
=================== ===================== =============== ================================
|
||||
|
||||
Important keywords
|
||||
++++++++++++++++++
|
||||
|
||||
|
||||
=================== ========================== ===================== =============== ================================
|
||||
keyword CONQUEST equivalence type default value description
|
||||
=================== ========================== ===================== =============== ================================
|
||||
``directory`` None ``str`` None directory used for storing input/output and calculation files
|
||||
``label`` None ``str`` None basename for working files (only used by ASE, eg. NEB)
|
||||
``kpts`` None ``list`` or ``tuple`` None k-points grid ; converted to CONQUEST Monkhorst-Pack grid
|
||||
``grid_cutoff`` ``Grid.GridCutoff`` ``float`` 100 integration grid in Ha
|
||||
``xc`` ``General.FunctionalType`` ``str`` 'PBE' exchange and correlation functional
|
||||
``self_consistent`` ``minE.SelfConsistent`` ``bool`` True choose either SCF or non-SCF
|
||||
``scf_tolerance`` ``minE.SCTolerance`` ``float`` 1e-6 Self-consistent-field convergence tolerance in Ha
|
||||
``nspin`` ``Spin.SpinPolarised`` ``int`` 1 spin polarisation: 1 for unpolarized or 2 for polarised
|
||||
``conquest_flags`` None ``dict`` None other CONQUET keyword arguments
|
||||
=================== ========================== ===================== =============== ================================
|
||||
|
||||
|
||||
Defaults keywords
|
||||
+++++++++++++++++
|
||||
|
||||
=============================== ========= =============== ================================
|
||||
keyword type default value description
|
||||
=============================== ========= =============== ================================
|
||||
``IO.WriteOutToASEFile`` ``bool`` True write ASE output file ; **must always be True when using ASE for post-processing**
|
||||
``IO.Iprint`` ``int`` 1 verbose for the output ; **must always be 1 when using ASE for post-processing**
|
||||
``DM.SolutionMethod`` ``str`` 'diagon' 'diagon' stands for diagonalisation other is 'ordern' (base on density matrix)
|
||||
``General.PseudopotentialType`` ``str`` 'Hamann' kind of pseudopotential other type are 'siesta' and 'abinit'
|
||||
``SC.MaxIters`` ``int`` 50 maximum number SCF cycles
|
||||
``AtomMove.TypeOfRun`` ``str`` 'static' 'static' stands for single (non)SCF other are 'md' or optimisation algorithms.
|
||||
``Diag.SmearingType`` ``int`` 1 1 for Methfessel-Paxton ; 0 for Fermi-Dirac
|
||||
``Diag.kT`` ``float`` 0.001 smearing temperature in Ha
|
||||
=============================== ========= =============== ================================
|
||||
|
||||
..
|
||||
``io.fractionalatomiccoords`` ``bool`` True atomic coordinates format for the structure file (fractional or cartesian)
|
||||
``basis.basisset`` ``str`` 'PAOs' type of basis set ; always 'PAOs' with ASE
|
||||
|
||||
Some examples
|
||||
-------------
|
||||
|
||||
An example of more advanced Calculator setup is given below for a SCF calculation on BCC-Na
|
||||
where for a PBE calculation using a k-point grid of :math:`6\times 6\times6` using the Fermi-Dirac
|
||||
distribution for the occupation with a smearing of 0.005 Ha::
|
||||
|
||||
struct = bulk('Na', crystalstructure='bcc', a=4.17, cubic=True)
|
||||
basis = {'Na' : {'file' : 'NaCQ.ion'}}
|
||||
|
||||
conquest_flags = {'Diag.SmearingType': 0,
|
||||
'Diag.kT' : 0.005}
|
||||
|
||||
struct.calc = Conquest(directory = 'Na_bcc_example',
|
||||
grid_cutoff = 90.0,
|
||||
self_consistent= True,
|
||||
xc = 'PBE',
|
||||
basis = basis,
|
||||
kpts = [6,6,6],
|
||||
nspin = 1,
|
||||
**conquest_flags)
|
||||
|
||||
struct.get_potential_energy()
|
||||
|
||||
..
|
||||
A dictionary containing a small number of mandatory keywords, listed below:
|
||||
|
||||
::
|
||||
|
||||
default_parameters = {
|
||||
'grid_cutoff' : 100, # DFT defaults
|
||||
'kpts' : None,
|
||||
'xc' : 'PBE',
|
||||
'scf_tolerance' : 1.0e-6,
|
||||
'nspin' : 1,
|
||||
'general.pseudopotentialtype' : 'Hamann', # CONQUEST defaults
|
||||
'basis.basisset' : 'PAOs',
|
||||
'io.iprint' : 2,
|
||||
'io.fractionalatomiccoords' : True,
|
||||
'mine.selfconsistent' : True,
|
||||
'sc.maxiters' : 50,
|
||||
'atommove.typeofrun' : 'static',
|
||||
'dm.solutionmethod' : 'diagon'}
|
||||
|
||||
The first five key/value pairs are special DFT parameters, the grid cutoff, the
|
||||
k-point mesh, the exchange-correlation functional, the SCF tolerance and the
|
||||
number of spins respectively. The rest are CONQUEST-specific input flags.
|
||||
|
||||
..
|
||||
The atomic species blocks are handled slightly differently, with a dictionary of
|
||||
their own. If the ``.ion`` files are present in the calculation directory, they
|
||||
can be specified as follows:
|
||||
|
||||
::
|
||||
|
||||
basis = {"H": {"valence_charge": 1.0,
|
||||
"number_of_supports": 1,
|
||||
"support_fn_range": 6.9},
|
||||
"O": {"valence_charge": 6.0,
|
||||
"number_of_supports": 4,
|
||||
"support_fn_range": 6.9}}
|
||||
|
||||
If the basis set ``.ion`` files are present in the directory containing the ASE
|
||||
script are pressent and are named ``element.ion``, then the relevant parameters
|
||||
will be parsed from the ``.ion`` files and included when the input file is
|
||||
written and this dictionary can be omitted. It is more important when, for
|
||||
example, setting up a multisite calculation, when the number of contracted
|
||||
support functions is different from the number in the ``.ion`` file.
|
||||
|
||||
ASE can also invoke the CONQUEST basis set generation tool, although care should
|
||||
be taken when generating basis sets:
|
||||
|
||||
::
|
||||
|
||||
basis = {"H": {"basis_size": "minimal",
|
||||
"pseudopotential_type": hamann",
|
||||
"gen_basis": True},
|
||||
"O": {"basis_size": "minimal",
|
||||
"pseudopotential_type": hamann",
|
||||
"gen_basis": True}}
|
||||
|
||||
Finally, defaults and other input flags can be defined in a new dictionary, and
|
||||
passed as an expanded set of keyword arguments.
|
||||
::
|
||||
|
||||
conquest_flags = {'DM.SolutionMethod' : 'ordern',
|
||||
'DM.L_range' : 12.0,
|
||||
'minE.LTolerance' : 1.0e-6}
|
||||
|
||||
Here is an example, combining the above. We set up a cubic diamond cell
|
||||
containing 8 atoms, and perform a single point energy calculation using the
|
||||
order(*N*) method (the default is diagonalisation, so we must specify all of the
|
||||
order(*N*) flags). We don't define a basis set, instead providing keywords that
|
||||
specify that a minimal basis set should be constructed using the MakeIonFiles
|
||||
basis generation tool.
|
||||
|
||||
::
|
||||
|
||||
import os
|
||||
from ase.build import bulk
|
||||
from ase.calculators.conquest import Conquest
|
||||
|
||||
CQ_ROOT = 'PATH_TO_CONQUEST_DIRECTORY'
|
||||
|
||||
os.environ['ASE_CONQUEST_COMMAND'] = 'mpirun -np 4 '+CQ_ROOT+'/bin/Conquest'
|
||||
os.environ["CQ_GEN_BASIS_CMD"] = CQ_ROOT+'/bin/MakeIonFiles"
|
||||
os.environ['CQ_PP_PATH'] = CQ_ROOT+'/pseudo-and-pao/'
|
||||
|
||||
diamond = bulk('C', 'diamond', a=3.6, cubic=True) # The atoms object
|
||||
conquest_flags = {'DM.SolutionMethod' : 'ordern', # Conquest keywords
|
||||
'DM.L_range' : 12.0,
|
||||
'minE.LTolerance' : 1.0e-6}
|
||||
|
||||
basis = {'C': {'basis_size' : 'minimal', # Generate a minimal basis
|
||||
'gen_basis' : True}
|
||||
|
||||
calc = Conquest(grid_cutoff = 80, # Set the calculator keywords
|
||||
xc = 'PBE',
|
||||
self_consistent=True,
|
||||
basis = basis,
|
||||
nspin = 1,
|
||||
**conquest_flags)
|
||||
|
||||
diamond.set_calculator(calc) # attach the calculator to the atoms object
|
||||
energy = diamond.get_potential_energy() # calculate the potential energy
|
||||
|
||||
Go to :ref:`top <ext-tools>`.
|
||||
|
||||
.. _et_ase_mssf:
|
||||
|
||||
Multisite support functions
|
||||
+++++++++++++++++++++++++++
|
||||
|
||||
Multisite support functions require a few additional keywords in the atomic
|
||||
species block, which can be specified as follows:
|
||||
|
||||
::
|
||||
|
||||
basis = {'C': {"basis_size": 'medium',
|
||||
"gen_basis": True,
|
||||
"pseudopotential_type": "hamann",
|
||||
"Atom.NumberofSupports": 4,
|
||||
"Atom.MultisiteRange": 7.0,
|
||||
"Atom.LFDRange": 7.0}}
|
||||
|
||||
Note that we are constructing a DZP basis set (size medium) with 13 primitive
|
||||
support functions using ``MakeIonFiles``, and contracting it to multisite basis
|
||||
of 4 support functions. The calculation requires a few more input flags, which
|
||||
are specified in the ``other_keywords`` dictionary:
|
||||
|
||||
::
|
||||
|
||||
other_keywords = {"Basis.MultisiteSF": True,
|
||||
"Multisite.LFD": True,
|
||||
"Multisite.LFD.Min.ThreshE": 1.0e-7,
|
||||
"Multisite.LFD.Min.ThreshD": 1.0e-7,
|
||||
"Multisite.LFD.Min.MaxIteration": 150,
|
||||
}
|
||||
|
||||
Go to :ref:`top <ext-tools>`.
|
||||
|
||||
.. _et_ase_load_dm:
|
||||
|
||||
Loading the K/L matrix
|
||||
++++++++++++++++++++++
|
||||
|
||||
Most calculation that involve incrementally moving atoms (molecular dynamics,
|
||||
geometry optimisation, equations of state, nudged elastic band etc.) can be made
|
||||
faster by using the K or L matrix from a previous calculation as the initial
|
||||
guess for a subsequent calculation in which that atoms have been moved slightly.
|
||||
This can be achieved by first performing a single point calculation to generate
|
||||
the first K/L matrix, then adding the following keywords to the calculator:
|
||||
|
||||
::
|
||||
|
||||
other_keywords = {"General.LoadL": True,
|
||||
"SC.MakeInitialChargeFromK": True}
|
||||
|
||||
These keywords respectively cause the K or L matrix to be loaded from file(s)
|
||||
``Kmatrix.i**.p*****``, and the initial charge density to be constructed from
|
||||
this matrix. In all subsequent calculations, the K or L matrix will be written
|
||||
at the end of the calculation and used as the initial guess for the subsequent
|
||||
ionic step.
|
||||
|
||||
Go to :ref:`top <ext-tools>`.
|
||||
|
||||
.. _et_eos:
|
||||
|
||||
Equation of state
|
||||
+++++++++++++++++
|
||||
|
||||
The following code computes the equation of state of diamond by doing single
|
||||
point calculations on a uniform grid of the ``a`` lattice parameter. It then
|
||||
interpolates the equation of state and uses ``matplotlib`` to generate a plot.
|
||||
|
||||
::
|
||||
|
||||
import scipy as sp
|
||||
from ase.build import bulk
|
||||
from ase.io.trajectory import Trajectory
|
||||
from ase.calculators.conquest import Conquest
|
||||
|
||||
|
||||
# Construct a unit cell
|
||||
diamond = bulk('C', 'diamond', a=3.6, cubic=True)
|
||||
|
||||
basis = {'C': {"basis_size": 'minimal',
|
||||
"gen_basis": True}}
|
||||
|
||||
calc = Conquest(grid_cutoff = 50,
|
||||
xc = "PBE",
|
||||
basis = basis,
|
||||
kpts = [4,4,4])
|
||||
|
||||
diamond.set_calculator(calc)
|
||||
|
||||
cell = diamond.get_cell()
|
||||
traj = Trajectory('diamond.traj', 'w') # save all results to trajectory
|
||||
|
||||
for x in sp.linspace(0.95, 1.05, 5): # grid for equation of state
|
||||
diamond.set_cell(cell*x, scale_atoms=True)
|
||||
diamond.get_potential_energy()
|
||||
traj.write(diamond)
|
||||
|
||||
from ase.io import read
|
||||
from ase.eos import EquationOfState
|
||||
|
||||
configs = read('diamond.traj@0:5')
|
||||
volumes = [diamond.get_volume() for diamond in configs]
|
||||
energies = [diamond.get_potential_energy() for diamond in configs]
|
||||
eos = EquationOfState(volumes, energies)
|
||||
v0, e0, B = eos.fit()
|
||||
|
||||
import matplotlib
|
||||
eos.plot('diamond-eos.pdf') # Plot the equation of state
|
||||
|
||||
Go to :ref:`top <ase-conquest>`.
|
||||
|
||||
|
||||
|
|
@ -188,248 +188,26 @@ output will look something like this:
|
|||
|
||||
Go to :ref:`top <ext-tools>`.
|
||||
|
||||
.. _et_ase:
|
||||
|
||||
Atomic Simulation Environment (ASE)
|
||||
-----------------------------------
|
||||
|
||||
ASE is a set of Python tools for setting up, manipulating, running, visualizing
|
||||
and analyzing atomistic simulations. ASE contains a CONQUEST interface, so that
|
||||
it can be used to calculate energies, forces and stresses for calculations that
|
||||
CONQUEST can't do (yet). Detailed instructions on how to install and invoke it
|
||||
can be found on its `website <https://wiki.fysik.dtu.dk/ase/>`_, but we provide
|
||||
some details and examples for the CONQUEST interface here.
|
||||
.. _et_ase_cq:
|
||||
|
||||
Note that the script will need to set environmental variables specifying the
|
||||
locations of the CONQUEST executable ``Conquest``, and if required, the basis
|
||||
set generation executable ``MakeIonFiles`` and pseudopotential database.
|
||||
`ASE <https://wiki.fysik.dtu.dk/ase>`_ is a set of
|
||||
Python tools for setting up, manipulating, running, visualizing and analyzing
|
||||
atomistic simulations. ASE contains a CONQUEST interface, also
|
||||
called *Calculator* so that it can be used to calculate ``energies``, ``forces``
|
||||
and ``stresses`` as inputs to other calculations such as `Phonon <https://wiki.fysik.dtu.dk/ase/ase/phonons.html#module-ase.phonons>`_
|
||||
or `NEB <https://wiki.fysik.dtu.dk/ase/ase/neb.html#module-ase.neb>`_ that
|
||||
are not implemented in CONQUEST. ASE is a versatil tool to manage CONQUEST
|
||||
calculations without pain either:
|
||||
|
||||
::
|
||||
* in a **direct** way where pre-processing, calculation and post-processing are managed on-the-fly by ASE,
|
||||
* or in an **indirect** way where the calculation step is performed outside the workflow, ie. on a supercomputer.
|
||||
|
||||
import os
|
||||
|
||||
# The command to run CONQUEST in parallel
|
||||
os.environ["ASE_CONQUEST_COMMAND"] = "mpirun -np 4 /path/to/Conquest_master"
|
||||
# Path to a database of pseudopotentials (for basis generation tool)
|
||||
os.environ["CQ_PP_PATH"] = "~/Conquest/PPDB/"
|
||||
# Path to the basis generation tool executable
|
||||
os.environ["CQ_GEN_BASIS_CMD"] = "/path/to/MakeIonFiles"
|
||||
The ASE repository containing the Conquest calculator can be
|
||||
found `here <https://gitlab.com/lionelalexandre/ase-Conquest/-/tree/master?ref_type=heads>`_.
|
||||
Detailed documentation on how to manage Conquest calculations
|
||||
with ASE is available :ref:`here <ase-conquest>`.
|
||||
|
||||
Go to :ref:`top <ext-tools>`.
|
||||
|
||||
.. _et_ase_input:
|
||||
|
||||
Keywords for generating the Conquest_input file
|
||||
+++++++++++++++++++++++++++++++++++++++++++++++
|
||||
|
||||
The calculator object contains a dictionray containing a small number of
|
||||
mandatory keywords, listed below:
|
||||
|
||||
::
|
||||
|
||||
default_parameters = {
|
||||
'grid_cutoff' : 100, # DFT defaults
|
||||
'kpts' : None,
|
||||
'xc' : 'PBE',
|
||||
'scf_tolerance' : 1.0e-6,
|
||||
'nspin' : 1,
|
||||
'general.pseudopotentialtype' : 'Hamann', # CONQUEST defaults
|
||||
'basis.basisset' : 'PAOs',
|
||||
'io.iprint' : 2,
|
||||
'io.fractionalatomiccoords' : True,
|
||||
'mine.selfconsistent' : True,
|
||||
'sc.maxiters' : 50,
|
||||
'atommove.typeofrun' : 'static',
|
||||
'dm.solutionmethod' : 'diagon'}
|
||||
|
||||
The first five key/value pairs are special DFT parameters, the grid cutoff, the
|
||||
k-point mesh, the exchange-correlation functional, the SCF tolerance and the
|
||||
number of spins respectively. The rest are CONQUEST-specific input flags.
|
||||
|
||||
The atomic species blocks are handled slightly differently, with a dictionary of
|
||||
their own. If the ``.ion`` files are present in the calculation directory, they
|
||||
can be specified as follows:
|
||||
|
||||
::
|
||||
|
||||
basis = {"H": {"valence_charge": 1.0,
|
||||
"number_of_supports": 1,
|
||||
"support_fn_range": 6.9},
|
||||
"O": {"valence_charge": 6.0,
|
||||
"number_of_supports": 4,
|
||||
"support_fn_range": 6.9}}
|
||||
|
||||
If the basis set ``.ion`` files are present in the directory containing the ASE
|
||||
script are pressent and are named ``element.ion``, then the relevant parameters
|
||||
will be parsed from the ``.ion`` files and included when the input file is
|
||||
written and this dictionary can be omitted. It is more important when, for
|
||||
example, setting up a multisite calculation, when the number of contracted
|
||||
support functions is different from the number in the ``.ion`` file.
|
||||
|
||||
ASE can also invoke the CONQUEST basis set generation tool, although care should
|
||||
be taken when generating basis sets:
|
||||
|
||||
::
|
||||
|
||||
basis = {"H": {"basis_size": "minimal",
|
||||
"pseudopotential_type": hamann",
|
||||
"gen_basis": True},
|
||||
"O": {"basis_size": "minimal",
|
||||
"pseudopotential_type": hamann",
|
||||
"gen_basis": True}}
|
||||
|
||||
Finally, non-mandatory input flags can be defined in a new dictionary, and
|
||||
passed as an expanded set of keyword arguments.
|
||||
|
||||
::
|
||||
|
||||
conquest_flags = {'IO.Iprint' : 1, # CONQUEST keywords
|
||||
'DM.SolutionMethod' : 'ordern',
|
||||
'DM.L_range' : 8.0,
|
||||
'minE.LTolerance' : 1.0e-6}
|
||||
|
||||
Here is an example, combining the above. We set up a cubic diamond cell
|
||||
containing 8 atoms, and perform a single point energy calculation using the
|
||||
order(N) method (the default is diagonalisation, so we must specify all of the
|
||||
order(N) flags). We don't define a basis set, instead providing keywords that
|
||||
specify that a minimal basis set should be constructed using the MakeIonFiles
|
||||
basis generation tool.
|
||||
|
||||
::
|
||||
|
||||
from ase.build import bulk
|
||||
from ase.calculators.conquest import Conquest
|
||||
|
||||
os.environ["ASE_CONQUEST_COMMAND"] = "mpirun -np 4 Conquest_master"
|
||||
os.environ["CQ_PP_PATH"] = "/Users/zamaan/Conquest/PPDB/"
|
||||
os.environ["CQ_GEN_BASIS_CMD"] = "MakeIonFiles"
|
||||
|
||||
diamond = bulk('C', 'diamond', a=3.6, cubic=True) # The atoms object
|
||||
conquest_flags = {'IO.Iprint' : 1, # Conquest keywords
|
||||
'DM.SolutionMethod' : 'ordern',
|
||||
'DM.L_range' : 8.0,
|
||||
'minE.LTolerance' : 1.0e-6}
|
||||
basis = {'C': {"basis_size" : 'minimal', # Generate a minimal basis
|
||||
"gen_basis" : True,
|
||||
"pseudopotential_type" : "hamann"}}
|
||||
|
||||
calc = Conquest(grid_cutoff = 80, # Set the calculator keywords
|
||||
xc="LDA",
|
||||
self_consistent=True,
|
||||
basis=basis,
|
||||
nspin=1,
|
||||
**conquest_flags)
|
||||
diamond.set_calculator(calc) # attach the calculator to the atoms object
|
||||
energy = diamond.get_potential_energy() # calculate the potential energy
|
||||
|
||||
Go to :ref:`top <ext-tools>`.
|
||||
|
||||
.. _et_ase_mssf:
|
||||
|
||||
Multisite support functions
|
||||
+++++++++++++++++++++++++++
|
||||
|
||||
Multisite support functions require a few additional keywords in the atomic
|
||||
species block, which can be specified as follows:
|
||||
|
||||
::
|
||||
|
||||
basis = {'C': {"basis_size": 'medium',
|
||||
"gen_basis": True,
|
||||
"pseudopotential_type": "hamann",
|
||||
"Atom.NumberofSupports": 4,
|
||||
"Atom.MultisiteRange": 7.0,
|
||||
"Atom.LFDRange": 7.0}}
|
||||
|
||||
Note that we are constructing a DZP basis set (size medium) with 13 primitive
|
||||
support functions using ``MakeIonFiles``, and contracting it to multisite basis
|
||||
of 4 support functions. The calculation requires a few more input flags, which
|
||||
are specified in the ``other_keywords`` dictionary:
|
||||
|
||||
::
|
||||
|
||||
other_keywords = {"Basis.MultisiteSF": True,
|
||||
"Multisite.LFD": True,
|
||||
"Multisite.LFD.Min.ThreshE": 1.0e-7,
|
||||
"Multisite.LFD.Min.ThreshD": 1.0e-7,
|
||||
"Multisite.LFD.Min.MaxIteration": 150,
|
||||
}
|
||||
|
||||
Go to :ref:`top <ext-tools>`.
|
||||
|
||||
.. _et_ase_load_dm:
|
||||
|
||||
Loading the K/L matrix
|
||||
++++++++++++++++++++++
|
||||
|
||||
Most calculation that involve incrementally moving atoms (molecular dynamics,
|
||||
geometry optimisation, equations of state, nudged elastic band etc.) can be made
|
||||
faster by using the K or L matrix from a previous calculation as the initial
|
||||
guess for a subsequent calculation in which that atoms have been moved slightly.
|
||||
This can be achieved by first performing a single point calculation to generate
|
||||
the first K/L matrix, then adding the following keywords to the calculator:
|
||||
|
||||
::
|
||||
|
||||
other_keywords = {"General.LoadL": True,
|
||||
"SC.MakeInitialChargeFromK": True}
|
||||
|
||||
These keywords respectively cause the K or L matrix to be loaded from file(s)
|
||||
``Kmatrix.i**.p*****``, and the initial charge density to be constructed from
|
||||
this matrix. In all subsequent calculations, the K or L matrix will be written
|
||||
at the end of the calculation and used as the initial guess for the subsequent
|
||||
ionic step.
|
||||
|
||||
Go to :ref:`top <ext-tools>`.
|
||||
|
||||
.. _et_eos:
|
||||
|
||||
Equation of state
|
||||
+++++++++++++++++
|
||||
|
||||
The following code computes the equation of state of diamond by doing single
|
||||
point calculations on a uniform grid of the ``a`` lattice parameter. It then
|
||||
interpolates the equation of state and uses ``matplotlib`` to generate a plot.
|
||||
|
||||
::
|
||||
|
||||
import scipy as sp
|
||||
from ase.build import bulk
|
||||
from ase.io.trajectory import Trajectory
|
||||
from ase.calculators.conquest import Conquest
|
||||
|
||||
|
||||
# Construct a unit cell
|
||||
diamond = bulk('C', 'diamond', a=3.6, cubic=True)
|
||||
|
||||
basis = {'C': {"basis_size": 'minimal',
|
||||
"gen_basis": True,
|
||||
"pseudopotential_type": "hamann"}}
|
||||
calc = Conquest(grid_cutoff = 50,
|
||||
xc = "LDA",
|
||||
basis = basis,
|
||||
kpts = [4,4,4]}
|
||||
diamond.set_calculator(calc)
|
||||
|
||||
cell = diamond.get_cell()
|
||||
traj = Trajectory('diamond.traj', 'w') # save all results to trajectory
|
||||
|
||||
for x in sp.linspace(0.95, 1.05, 5): # grid for equation of state
|
||||
diamond.set_cell(cell*x, scale_atoms=True)
|
||||
diamond.get_potential_energy()
|
||||
traj.write(diamond)
|
||||
|
||||
from ase.io import read
|
||||
from ase.eos import EquationOfState
|
||||
|
||||
configs = read('diamond.traj@0:5')
|
||||
volumes = [diamond.get_volume() for diamond in configs]
|
||||
energies = [diamond.get_potential_energy() for diamond in configs]
|
||||
eos = EquationOfState(volumes, energies)
|
||||
v0, e0, B = eos.fit()
|
||||
|
||||
import matplotlib
|
||||
eos.plot('diamond-eos.pdf') # Plot the equation of state
|
||||
|
||||
Go to :ref:`top <ext-tools>`.
|
||||
|
||||
|
|
|
|||
|
|
@ -49,6 +49,7 @@ User Guide
|
|||
* :doc:`moldyn`
|
||||
* :doc:`post-proc`
|
||||
* :doc:`ext-tools`
|
||||
* :doc:`ase-conquest`
|
||||
* :doc:`errors`
|
||||
* :doc:`input_tags`
|
||||
|
||||
|
|
@ -65,6 +66,7 @@ User Guide
|
|||
strucrelax
|
||||
moldyn
|
||||
post-proc
|
||||
ase-conquest
|
||||
ext-tools
|
||||
errors
|
||||
input_tags
|
||||
|
|
@ -93,6 +95,7 @@ Tutorials
|
|||
Theory
|
||||
------
|
||||
|
||||
* :doc:`theory_energy_force_stress`
|
||||
* :doc:`theory-strucrelax`
|
||||
* :doc:`theory-md`
|
||||
|
||||
|
|
@ -101,6 +104,7 @@ Theory
|
|||
:hidden:
|
||||
:caption: Theory
|
||||
|
||||
theory_energy_force_stress
|
||||
theory-strucrelax
|
||||
theory-md
|
||||
|
||||
|
|
|
|||
|
|
@ -1403,6 +1403,16 @@ General.GapThreshold (*real*)
|
|||
General.only_Dispersion (*boolean*)
|
||||
Selects only DFT\_D2 calculation (no electronic structure etc)
|
||||
|
||||
General.MixXCGGAInOut (*real*)
|
||||
For non-SCF calculations only, chooses how to mix the proportions of
|
||||
GGA XC stress contribution (from the change of the electron density
|
||||
gradient) found using input (0.0 gives pure input) and output (1.0
|
||||
gives pure output) densities. Note that this is an approximation but
|
||||
varying the value significantly away from 0.5 will give inconsistency
|
||||
between stress and energy.
|
||||
|
||||
*default*: 0.5
|
||||
|
||||
Go to :ref:`top <input_tags>`.
|
||||
|
||||
.. _advanced_atomic_spec_tags:
|
||||
|
|
|
|||
|
|
@ -211,6 +211,8 @@ in a well-converged input charge density: set ``minE.SelfConsistent
|
|||
F`` and ``General.LoadRho T`` (which will require that the converged
|
||||
charge density is written out by CONQUEST by setting ``IO.DumpChargeDensity T``).
|
||||
|
||||
Go to :ref:`top <post-proc>`.
|
||||
|
||||
.. _pp_pDOS:
|
||||
|
||||
Atom-projected DOS
|
||||
|
|
@ -293,5 +295,62 @@ also be specified:
|
|||
where the final tag sets the minimum and maximum values relative to
|
||||
the Fermi level.
|
||||
|
||||
If you only want to produce pDOS for a few atoms, then you can set
|
||||
the variable ``Process.n_atoms_pDOS`` and list the atoms you want
|
||||
in the block ``pDOS_atoms``:
|
||||
|
||||
::
|
||||
|
||||
Process.n_atoms_pDOS 2
|
||||
%block pDOS_atoms
|
||||
1
|
||||
12
|
||||
%endblock
|
||||
|
||||
Go to :ref:`top <post-proc>`.
|
||||
|
||||
.. _pp_band_str
|
||||
|
||||
Band structure
|
||||
++++++++++++++
|
||||
|
||||
The band structure of a material can be generated by CONQUEST by performing
|
||||
a non-self-consistent calculation, after reading a well-converged charge density:
|
||||
set ``minE.SelfConsistent F`` and ``General.LoadRho T`` (remember that to write
|
||||
a converged charge density from CONQUEST you set ``IO.DumpChargeDensity T``).
|
||||
The k-points required can be specified as lines of points in k-space;
|
||||
setting ``Diag.KspaceLines T`` enables this (replacing the usual MP mesh), while the number of lines
|
||||
(e.g. Gamma to L; L to X; would be two lines) is set with ``Diag.NumKptLines``
|
||||
and the number of points along a line with ``Diag.NumKpts``. The k-point lines
|
||||
themselves are set with a block labelled ``Diag.KpointLines`` which should have
|
||||
two entries (starting and finishing k-points) for each k-point line.
|
||||
(In constructing the k-point list, CONQUEST will automatically remove any duplicate
|
||||
points, so that the output can be plotted smoothly.) So to create
|
||||
a bandstructure from X-:math:`\Gamma`-L-X (3 lines: X-:math:`\Gamma`; :math:`\Gamma`-L;
|
||||
L-X) with 11 points in each line, you would use the following input:
|
||||
|
||||
::
|
||||
|
||||
Diag.KspaceLines T
|
||||
Diag.NumKptLines 3
|
||||
Diag.NumKpts 11
|
||||
%block Diag.KpointLines
|
||||
0.5 0.0 0.0
|
||||
0.0 0.0 0.0
|
||||
0.0 0.0 0.0
|
||||
0.5 0.5 0.5
|
||||
0.5 0.5 0.5
|
||||
0.5 0.0 0.0
|
||||
%endblock
|
||||
|
||||
After running CONQUEST, setting ``Process.Job bst`` and running the post-processing
|
||||
will read the resulting
|
||||
``eigenvalues.dat`` file, and produce a file ``BandStructure.dat``. The x-axis will
|
||||
be the k-point index by default, but specifying ``Process.BandStrucAxis`` (taking
|
||||
value ``n`` for index, ``x``, ``y`` or ``z`` for a single direction in k-space,
|
||||
or ``a`` to give all k-point coordinates) will allow you to control this. Limits
|
||||
on the energies to select the bands produced can be set
|
||||
with ``Process.min_DOS_E`` and ``Process.max_DOS_E``.
|
||||
|
||||
Go to :ref:`top <post-proc>`.
|
||||
|
||||
|
|
|
|||
|
|
@ -7,6 +7,30 @@
|
|||
|
||||
%% Saved with string encoding Unicode (UTF-8)
|
||||
|
||||
@article{Torralba:2009nr,
|
||||
author = {Torralba, Antonio S and Bowler, David R and Miyazaki, Tsuyoshi and Gillan, Michael J},
|
||||
doi = {10.1021/ct8005425},
|
||||
journal = {J. Chem. Theory Comput.},
|
||||
pages = {1499},
|
||||
title = {{Non-self-consistent Density-Functional Theory Exchange-Correlation Forces for GGA Functionals}},
|
||||
volume = {5},
|
||||
year = {2009}
|
||||
}
|
||||
|
||||
@article{Louie:1982aa,
|
||||
title = {Nonlinear ionic pseudopotentials in spin-density-functional calculations},
|
||||
author = {Louie, Steven G. and Froyen, Sverre and Cohen, Marvin L.},
|
||||
journal = {Phys. Rev. B},
|
||||
volume = {26},
|
||||
issue = {4},
|
||||
pages = {1738--1742},
|
||||
numpages = {0},
|
||||
year = {1982},
|
||||
month = {Aug},
|
||||
publisher = {American Physical Society},
|
||||
doi = {10.1103/PhysRevB.26.1738}
|
||||
}
|
||||
|
||||
@article{Resta:1992aa,
|
||||
Author = {R. Resta},
|
||||
Doi = {10.1002/(SICI)1097-461X(1999)75:4/5%3C599::AID-QUA25%3E3.0.CO;2-8},
|
||||
|
|
|
|||
|
|
@ -0,0 +1,183 @@
|
|||
.. _theory_energy_force_stress:
|
||||
|
||||
=======================================
|
||||
Background on energy, forces and stress
|
||||
=======================================
|
||||
|
||||
A number of different ways of formulating the energy exist in Conquest at the moment, involving both self-consistent and non-self-consistent densities and potentials, both with and without the neutral atom potential. With self-consistency, all formulations should give the same result, though numerical issues may give small differences; without self-consistency the Harris-Foulkes functional is more accurate.
|
||||
|
||||
Self-consistent calculations
|
||||
~~~~~~~~~~~~~~~~~~~~~~~~~~~~
|
||||
|
||||
We define the energy in Conquest in two ways that are equivalent at the self-consistent ground state. The Harris-Foulkes energy is given as:
|
||||
|
||||
.. math::
|
||||
E_{HF} = 2\mathrm{Tr}\left[KH\right] + \Delta E_{Ha} + \Delta E_{XC} + E_{II}
|
||||
|
||||
where the first term is the band structure energy, equivalent to the sum over the energies of the occupied states, the second two terms compensate for double counting and the final term gives the ion-ion interaction:
|
||||
|
||||
.. math::
|
||||
E_{II} = \frac{1}{2}\left( \sum_{ij} \frac{Z_i Z_j}{\mid \mathbf{R}_i - \mathbf{R}_j \mid} \right)
|
||||
|
||||
The Hamiltonian is defined as:
|
||||
|
||||
.. math::
|
||||
\hat{H} = \hat{T} + \hat{V}_{L} + \hat{V}_{NL} + V_{Ha} + V_{XC}
|
||||
|
||||
where the operators are the kinetic energy, the local and non-local pseudopotentials, the Hartree potential, defined as :math:`V_{Ha} = \int d\mathbf{r}^\prime n(\mathbf{r}^\prime)/\mid \mathbf{r} - \mathbf{r}^\prime\mid`, and the exchange-correlation potential. The alternative form, often known as the DFT energy, is:
|
||||
|
||||
.. math::
|
||||
E_{DFT} = 2\mathrm{Tr}\left[K(T + V_{L} + V_{NL})\right] + E_{Ha} + E_{XC} + E_{II}
|
||||
|
||||
with the Hartree energy defined as usual:
|
||||
|
||||
.. math::
|
||||
E_{Ha} = \frac{1}{2}\int\int d\mathbf{r}d\mathbf{r}^{\prime} \frac{n(\mathbf{r})n(\mathbf{r}^\prime)}{\mid \mathbf{r} - \mathbf{r}^\prime\mid}
|
||||
|
||||
along with the exchange-correlation energy:
|
||||
|
||||
.. math::
|
||||
E_{XC} = \int d\mathbf{r} \epsilon_{XC}\left[n\right] n(\mathbf{r})
|
||||
|
||||
For the Harris-Foulkes and DFT energies to be equal, it is easy to see that the double counting correction terms in the Harris-Foulkes formalism must be:
|
||||
|
||||
.. math::
|
||||
\Delta E_{Ha} = -E_{Ha} = -\frac{1}{2}\int\int d\mathbf{r}d\mathbf{r}^{\prime} \frac{n(\mathbf{r})n(\mathbf{r}^\prime)}{\mid \mathbf{r} - \mathbf{r}^\prime\mid}
|
||||
|
||||
and
|
||||
|
||||
.. math::
|
||||
\Delta E_{XC} = \int d\mathbf{r} \left(\epsilon_{XC}[n] - V_{XC}[n]\right)n(\mathbf{r})
|
||||
|
||||
When calculating forces and stress with self-consistency, we generally use the differentials of the DFT energy rather than the Harris-Foulkes energy; this enables us to separate contributions that are calculated in different ways (in particular on those that are calculated on the integration grid from those that are not).
|
||||
|
||||
Go to :ref:`top <theory_energy_force_stress>`.
|
||||
|
||||
.. _th_efs_nap:
|
||||
|
||||
Neutral atom potential
|
||||
----------------------
|
||||
|
||||
In a DFT code using local orbitals as basis functions, the total energy is most conveniently written in terms of the interaction of neutral atoms: this is simply a reformulation of the total energy which, in particular, reduces the ion-ion interaction to a sum over short-range pair-wise interactions. The charge density of interest is now the *difference* between the total charge density and a superposition of atomic densities, notated as :math:`\delta n(\mathbf{r}) = n(\mathbf{r}) - \sum_i n_i(\mathbf{r})` for atomic densities :math:`n_i(\mathbf{r})`. We write:
|
||||
|
||||
.. math::
|
||||
E_{DFT, NA} = 2\mathrm{Tr}\left[K(T + V_{NA} + V_{NL})\right] + E_{\delta Ha} + E_{XC} + E_{SII}
|
||||
|
||||
where the second term is defined as:
|
||||
|
||||
.. math::
|
||||
E_{\delta Ha} = \frac{1}{2}\int\int d\mathbf{r}d\mathbf{r}^{\prime} \frac{\delta n(\mathbf{r})\delta n(\mathbf{r}^\prime)}{\mid \mathbf{r} - \mathbf{r}^\prime\mid}
|
||||
|
||||
The final term, the screened ion-ion interaction, is short-ranged, and written as:
|
||||
|
||||
.. math::
|
||||
E_{SII} = \frac{1}{2}\left( \sum_{ij} \frac{Z_i Z_j}{\mid \mathbf{R}_i - \mathbf{R}_j \mid} - \int d\mathbf{r} n_i(\mathbf{r})V_{Ha,j}(\mathbf{r}) \right)
|
||||
|
||||
where :math:`V_{Ha,i}(\mathbf{r})` is the Hartree potential from the atomic density :math:`n_i(\mathbf{r})`. We define the neutral atom potential for an atom as :math:`V_{NA,i}(\mathbf{r}) = V_{L,i}(\mathbf{r}) + V_{Ha,i}(\mathbf{r})`, combining the local potential and the Hartree potential for the atomic density; the overall neutral atom potential is given as the sum over the atomic densities, :math:`V_{NA}(\mathbf{r}) = \sum_i V_{NA,i}(\mathbf{r})`. If we write the pseudo-atomic density as :math:`n_{PAD}(\mathbf{r}) = \sum_i n_i(\mathbf{r})` then we can also write :math:`V_{NA}(\mathbf{r}) = V_L(\mathbf{r}) + V_{Ha, PAD}(\mathbf{r})`.
|
||||
|
||||
In this case, we can write the Harris-Foulkes energy as:
|
||||
|
||||
.. math::
|
||||
E_{HF} = 2\mathrm{Tr}\left[KH\right] + \Delta E_{Ha} + \Delta E_{XC} + E_{SII}
|
||||
|
||||
with the Hamiltonian defined as:
|
||||
|
||||
.. math::
|
||||
\hat{H} = \hat{T} + \hat{V}_{NA} + \hat{V}_{NL} + V_{\delta Ha} + V_{XC}
|
||||
|
||||
where :math:`V_{\delta Ha}(\mathbf{r}) = \int d\mathbf{r^\prime} \delta n(\mathbf{r^\prime})/\mid \mathbf{r} - \mathbf{r}^\prime\mid`. Accordingly, the double counting Hartree correction term has to change:
|
||||
|
||||
.. math::
|
||||
\Delta E_{Ha} = -E_{\delta Ha} - \int d\mathbf{r} \delta n(\mathbf{r})\sum_i V_{Ha,i}(\mathbf{r}).
|
||||
|
||||
Go to :ref:`top <theory_energy_force_stress>`.
|
||||
|
||||
.. _th_efs_nsc:
|
||||
|
||||
Non-self-consistent calculations
|
||||
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
|
||||
|
||||
In non-self-consistent calculations, we use the Harris-Foulkes functional, along with a reasonable guess for the input density, which is normally taken as the superposition of atomic densities, :math:`n_{in}(\mathbf{r})` and write:
|
||||
|
||||
.. math::
|
||||
E_{NSC} = 2\mathrm{Tr}\left[KH\right] + \Delta E_{Ha}\left[n_{in}\right] + \Delta E_{XC}\left[n_{in}\right] + E_{II}
|
||||
|
||||
Notice that we effectively have two densities being used here: :math:`n_{in}` (which is normally the superposition of atomic densities used in the neutral atom case) and and *effective* output density, :math:`n_{out} = \sum_{ij} \phi_i K_{ij} \phi_j` which comes from the band energy (first term); this complicates the calculation of forces and stress compared to the self-consistent case, as we have to consider contributions from both densities.
|
||||
|
||||
For the :ref:`neutral atom potential <th_efs_nap>`, :math:`\delta n(\mathbf{r}) = 0` by definition, which also means that :math:`E_{\delta Ha} = 0` and :math:`\Delta E_{Ha}=0`.
|
||||
|
||||
Go to :ref:`top <theory_energy_force_stress>`.
|
||||
|
||||
.. _th_efs_pcc:
|
||||
|
||||
Partial core corrections
|
||||
~~~~~~~~~~~~~~~~~~~~~~~~
|
||||
|
||||
Also known as non-linear core corrections, partial core corrections (PCC) :cite:`efs-Louie:1982aa` add a model core charge to the pseudopotential to allow for the non-linear exchange-correlation interation between core and valence charge (which is linearised in standard pseudopotentials); this generally improves the accuracy of the pseudopotential. The exchange-correlation potential is evaluated in terms of the combined charge density, :math:`n_v(\mathbf{r}) + n_c(\mathbf{r})` where the valence charge is input or output charge density defined above: :math:`V_{XC}\left[ n_v + n_c \right]`. The exchange-correlation energy becomes:
|
||||
|
||||
.. math::
|
||||
E_{XC} = \int d\mathbf{r} \left(n_v(\mathbf{r}) + n_c(\mathbf{r})\right) V_{XC}\left[ n_v + n_c \right] .
|
||||
|
||||
Once this change to the charge density has been made, there is no change to the DFT energy. However, the double counting term for Harris-Foulkes needs redefining, since XC contribution to the band energy is :math:`2Tr[KV_{XC}] = \int d\mathbf{r} n_v(\mathbf{r}) V_{XC}[n_v + n_c]`. We write:
|
||||
|
||||
.. math::
|
||||
\Delta E_{XC} &=& \int d\mathbf{r} \left(n_v(\mathbf{r}) + n_c(\mathbf{r})\right)\epsilon_{XC}[n_v + n_c] - \int d\mathbf{r} n_v(\mathbf{r})V_{XC}[n_v + n_c]\\
|
||||
&=& \int d\mathbf{r} n_c(\mathbf{r})\epsilon_{XC}[n_v + n_c] + \int d\mathbf{r} \left(\epsilon_{XC}[n_v + n_c] - V_{XC}[n_v + n_c]\right)n_v(\mathbf{r})
|
||||
|
||||
There is an extra factor of :math:`\int d\mathbf{r} n_c(\mathbf{r})\epsilon_{XC}[n_v + n_c]` over and above the usual term.
|
||||
|
||||
Go to :ref:`top <theory_energy_force_stress>`.
|
||||
|
||||
.. _th_efs_fs:
|
||||
|
||||
Forces and Stresses
|
||||
~~~~~~~~~~~~~~~~~~~
|
||||
|
||||
It is important that the forces and stresses be the exact derivatives of the energy, for consistency. In particular, this means that as the energy is calculated in different ways for different contributions, the force or stress contribution must be calculated in the same way.
|
||||
|
||||
Go to :ref:`top <theory_energy_force_stress>`.
|
||||
|
||||
.. _th_efs_for:
|
||||
|
||||
Forces
|
||||
------
|
||||
|
||||
Forces are defined as the change in energy with respect to atomic positions; as the basis functions move with the atoms, these changes will also include Pulay terms. The forces found in Conquest are documented extensively elsewhere :cite:`efs-Miyazaki2004,efs-Torralba:2009nr` though the changes needed to account for :ref:`PCC <th_efs_pcc>`, particularly in the :ref:`non-self-consistent <th_efs_nsc>` case, have not been published and are given here for completeness. As well as the Hellmann-Feynman forces (which come from the movement of the local and non-local pseudopotentials with the atoms) we define Pulay forces (divided into two parts, known as phi-Pulay which come from changes in the Hamiltonian matrix, and S-Pulay, which come from changes in the overlap matrix; the phi-Pulay forces are calculated in three contributions, which depend on how the respective parts of the Hamiltonian matrix are calculated: the kinetic energy; the non-local pseudopotential; and the remaining terms which are all found on the integration grid). The ion-ion interactions also contribute forces.
|
||||
|
||||
The inclusion of :ref:`PCC <th_efs_pcc>` adds an extra term to the forces in all calculations, which comes from the change of the core density as the atoms move; the force on atom :math:`i` is given as:
|
||||
|
||||
.. math::
|
||||
\mathbf{F}^{PCC}_i = -\int d\mathbf{r} \nabla_i n^c_i(\mathbf{r}) V_{XC}[n_v + n_c]
|
||||
|
||||
If the :ref:`non-self-consistent <th_efs_nsc>` formalism is used, then a further term is added (the non-self-consistent force changes) to include the gradient of the core charge. The non-self-consistent force is now written as:
|
||||
|
||||
.. math::
|
||||
\mathbf{F}^{NSC}_i = -\int d\mathbf{r} V_{\delta Ha}(\mathbf{r}) \nabla_i n^v_i(\mathbf{r}) - \int d\mathbf{r} \delta n(\mathbf{r}) V_{XC}^\prime\left[n^{in}_v + n_c\right] \left( \nabla_i n^v_i(\mathbf{r}) + \nabla_i n^c_i(\mathbf{r}) \right)
|
||||
|
||||
where :math:`V^\prime_{XC}` is the derivative of the exchange-correlation potential with respect to charge density.
|
||||
|
||||
Go to :ref:`top <theory_energy_force_stress>`.
|
||||
|
||||
.. _th_efs_str:
|
||||
|
||||
Stress
|
||||
------
|
||||
|
||||
The stress includes all contributions to the change of energy with the lattice constants; the calculation of stress in Conquest is documented in a paper being prepared for publication, but we give a brief overview here. As Conquest uses orthorhombic cells, only the diagonal stress components (:math:`\sigma_{\alpha\alpha}`) are calculated.
|
||||
|
||||
In most cases, forces also contribute to the stress; it is easy to show that the stress contribution is given by:
|
||||
|
||||
.. math::
|
||||
\sigma_{\alpha\alpha} = \sum_i F_{i\alpha}R_{i\alpha}
|
||||
|
||||
where :math:`R_{i\alpha}` is the position of the atom. As well as these contributions, there are more subtle terms. Any energies calculated on the grid will contribute to the stress as the integration grid changes with cell size (the stress is simply the energy calculated), and the Hartree potential contributes a term related to the change in the reciprocal lattice vectors (as it is calculated by Fourier transforming the charge density). If the exchange-correlation functional is a GGA functional, then a further term coming from the change of the gradient of the density with the cell size arises. (For non-self-consistent calculations this leads to some complications, as this term technically requires both input and output densities; at present, we approximate this as a mixture of the term calculated with input density and the term calculated with output density; the proportion can be adjusted using the parameter ``General.MixXCGGAInOut`` documented in the :ref:`Advanced and obscure tags <advanced_general_tags>` section of the manual, though we do not recommend changing it.)
|
||||
|
||||
Go to :ref:`top <theory_energy_force_stress>`.
|
||||
|
||||
.. bibliography:: references.bib
|
||||
:cited:
|
||||
:labelprefix: EFS
|
||||
:keyprefix: efs-
|
||||
:style: unsrt
|
||||
|
||||
Go to :ref:`top <theory_energy_force_stress>`.
|
||||
|
|
@ -1235,8 +1235,8 @@ contains
|
|||
end if
|
||||
|
||||
! dump the L matrix if required
|
||||
if (n_dumpL>0 .and. mod (n_iter, n_dumpL) == 0) then
|
||||
call dump_pos_and_matrices(index=unit_DM_save)
|
||||
if (n_dumpL>0) then
|
||||
if(mod (n_iter, n_dumpL) == 0) call dump_pos_and_matrices(index=unit_DM_save)
|
||||
endif
|
||||
!if (flag_dump_L) then
|
||||
! if (mod (n_iter, n_dumpL) == 0) then
|
||||
|
|
|
|||
|
|
@ -140,6 +140,8 @@ contains
|
|||
real(double) :: rcut
|
||||
real(double) :: r1, r2, r3, r4, core_charge, gauss_charge
|
||||
real(double) :: val, theta, phi, r_tmp
|
||||
integer :: max_num_blocks, current_num_blocks
|
||||
real(double), dimension(:), allocatable :: temp_block_storage
|
||||
|
||||
call start_timer(tmr_std_basis)
|
||||
call start_timer(tmr_std_allocation)
|
||||
|
|
@ -157,6 +159,9 @@ contains
|
|||
no_of_ib_ia = 0
|
||||
gridfunctions(pao_fns)%griddata = zero
|
||||
|
||||
max_num_blocks = maxval(npao_species)*n_pts_in_block
|
||||
allocate(temp_block_storage(max_num_blocks))
|
||||
|
||||
! loop arround grid points in the domain, and for each
|
||||
! point, get the PAO values
|
||||
|
||||
|
|
@ -168,9 +173,6 @@ contains
|
|||
iatom=0
|
||||
do ipart=1,naba_atoms_of_blocks(atomf)%no_of_part(iblock)
|
||||
jpart=naba_atoms_of_blocks(atomf)%list_part(ipart,iblock)
|
||||
if(jpart > DCS_parts%mx_gcover) then
|
||||
call cq_abort('single_PAO_to_grid: JPART ERROR ',ipart,jpart)
|
||||
endif
|
||||
ind_part=DCS_parts%lab_cell(jpart)
|
||||
do ia=1,naba_atoms_of_blocks(atomf)%no_atom_on_part(ipart,iblock)
|
||||
iatom=iatom+1
|
||||
|
|
@ -178,15 +180,6 @@ contains
|
|||
icover= DCS_parts%icover_ibeg(jpart)+ii-1
|
||||
ig_atom= id_glob(parts%icell_beg(ind_part)+ii-1)
|
||||
|
||||
if(parts%icell_beg(ind_part) + ii-1 > ni_in_cell) then
|
||||
call cq_abort('single_PAO_to_grid: globID ERROR ', &
|
||||
ii,parts%icell_beg(ind_part))
|
||||
endif
|
||||
if(icover > DCS_parts%mx_mcover) then
|
||||
call cq_abort('single_PAO_to_grid: icover ERROR ', &
|
||||
icover,DCS_parts%mx_mcover)
|
||||
endif
|
||||
|
||||
xatom=DCS_parts%xcover(icover)
|
||||
yatom=DCS_parts%ycover(icover)
|
||||
zatom=DCS_parts%zcover(icover)
|
||||
|
|
@ -197,17 +190,21 @@ contains
|
|||
rcut = r_h + RD_ERR
|
||||
call check_block (xblock,yblock,zblock,xatom,yatom,zatom, rcut, & ! in
|
||||
npoint,ip_store,r_store,x_store,y_store,z_store,n_pts_in_block) !out
|
||||
r_from_i = sqrt((xatom-xblock)**2+(yatom-yblock)**2+ &
|
||||
(zatom-zblock)**2 )
|
||||
|
||||
if(npoint > 0) then
|
||||
!offset_position = (no_of_ib_ia-1) * npao * n_pts_in_block
|
||||
current_num_blocks = npao_species(the_species)*n_pts_in_block
|
||||
temp_block_storage = zero
|
||||
offset_position = no_of_ib_ia
|
||||
!$omp parallel do default(none) &
|
||||
!$omp schedule(dynamic) &
|
||||
!$omp reduction(+: temp_block_storage) &
|
||||
!$omp shared(n_pts_in_block, pao, the_species, offset_position, &
|
||||
!$omp npoint, no_of_ib_ia, ip_store, r_store, x_store, y_store, z_store) &
|
||||
!$omp private(r_from_i, ip, position, ipoint, &
|
||||
!$omp x, y, z, l1, acz, m1, count1, val)
|
||||
do ip=1,npoint
|
||||
ipoint=ip_store(ip)
|
||||
position= offset_position + ipoint
|
||||
if(position > gridfunctions(pao_fns)%size) call cq_abort &
|
||||
('single_pao_to_grid: position error ', position, gridfunctions(pao_fns)%size)
|
||||
|
||||
r_from_i = r_store(ip)
|
||||
x = x_store(ip)
|
||||
|
|
@ -219,21 +216,23 @@ contains
|
|||
do acz = 1,pao(the_species)%angmom(l1)%n_zeta_in_angmom
|
||||
do m1=-l1,l1
|
||||
call evaluate_pao(the_species,l1,acz,m1,x,y,z,val)
|
||||
if(position+(count1-1)*n_pts_in_block > gridfunctions(pao_fns)%size) &
|
||||
call cq_abort('single_pao_to_grid: position error ', &
|
||||
position, gridfunctions(pao_fns)%size)
|
||||
gridfunctions(pao_fns)%griddata(position+(count1-1)*n_pts_in_block) = val
|
||||
!gridfunctions(pao_fns)%griddata(position+(count1-1)*n_pts_in_block) = val
|
||||
temp_block_storage(ipoint + (count1-1)*n_pts_in_block) = val
|
||||
count1 = count1+1
|
||||
end do ! m1
|
||||
end do ! acz
|
||||
end do ! l1
|
||||
enddo ! ip=1,npoint
|
||||
!$omp end parallel do
|
||||
gridfunctions(pao_fns)%griddata(no_of_ib_ia+1:no_of_ib_ia+current_num_blocks) = &
|
||||
temp_block_storage(1:current_num_blocks)
|
||||
endif! (npoint > 0) then
|
||||
no_of_ib_ia = no_of_ib_ia + npao_species(the_species)*n_pts_in_block
|
||||
enddo ! naba_atoms
|
||||
enddo ! naba_part
|
||||
endif !(naba_atoms_of_blocks(atomf)%no_of_part(iblock) > 0) !naba atoms?
|
||||
enddo ! iblock : primary set of blocks
|
||||
deallocate(temp_block_storage)
|
||||
call my_barrier()
|
||||
call start_timer(tmr_std_allocation)
|
||||
deallocate(ip_store,x_store,y_store,z_store,r_store)
|
||||
|
|
@ -254,6 +253,7 @@ contains
|
|||
!! Used for gradients of energy wrt atomic coordinates
|
||||
!!
|
||||
!! This subroutine is based on sub:single_PAO_to_grid in PAO_grid_transform_module.f90.
|
||||
!! TODO: There is a lot of code duplication between single_PAO_to_grid and single_PAO_to_grad
|
||||
!!
|
||||
!! INPUTS
|
||||
!!
|
||||
|
|
@ -313,6 +313,8 @@ contains
|
|||
real(double) :: rcut
|
||||
real(double) :: r1, r2, r3, r4, core_charge, gauss_charge
|
||||
real(double) :: val, theta, phi, r_tmp
|
||||
integer :: max_num_blocks, current_num_blocks
|
||||
real(double), dimension(:), allocatable :: temp_block_storage
|
||||
|
||||
call start_timer(tmr_std_basis)
|
||||
call start_timer(tmr_std_allocation)
|
||||
|
|
@ -330,21 +332,12 @@ contains
|
|||
no_of_ib_ia = 0
|
||||
gridfunctions(pao_fns)%griddata = zero
|
||||
|
||||
max_num_blocks = maxval(npao_species)*n_pts_in_block
|
||||
allocate(temp_block_storage(max_num_blocks))
|
||||
|
||||
! loop arround grid points in the domain, and for each
|
||||
! point, get the d_PAO/d_R values
|
||||
|
||||
!$omp parallel do default(none) &
|
||||
!$omp schedule(dynamic) &
|
||||
!$omp shared(domain, dcellx_block, dcelly_block, dcellz_block, atomf, &
|
||||
!$omp naba_atoms_of_blocks, parts, r_h, n_pts_in_block, direction, &
|
||||
!$omp pao, pao_fns, npao_species, dcs_parts, id_glob, species_glob, &
|
||||
!$omp gridfunctions) &
|
||||
!$omp private(iblock, iatom, ipart, jpart, ind_part, xblock, yblock, zblock, &
|
||||
!$omp ii, ia, icover, ig_atom, xatom, yatom, zatom, the_species, &
|
||||
!$omp rcut, r_from_i, ip, position, offset_position, ipoint, &
|
||||
!$omp npoint, ip_store, r_store, x_store, y_store, z_store, x, y, z, &
|
||||
!$omp l1, acz, m1, count1, val) &
|
||||
!$omp firstprivate(no_of_ib_ia)
|
||||
do iblock = 1, domain%groups_on_node ! primary set of blocks
|
||||
xblock=(domain%idisp_primx(iblock)+domain%nx_origin-1)*dcellx_block
|
||||
yblock=(domain%idisp_primy(iblock)+domain%ny_origin-1)*dcelly_block
|
||||
|
|
@ -370,11 +363,19 @@ contains
|
|||
rcut = r_h + RD_ERR
|
||||
call check_block (xblock,yblock,zblock,xatom,yatom,zatom, rcut, & ! in
|
||||
npoint,ip_store,r_store,x_store,y_store,z_store,n_pts_in_block) !out
|
||||
r_from_i = sqrt( (xatom-xblock)**2 + (yatom-yblock)**2 + (zatom-zblock)**2 )
|
||||
|
||||
if(npoint > 0) then
|
||||
!offset_position = (no_of_ib_ia-1) * npao * n_pts_in_block
|
||||
! Temporary storage
|
||||
current_num_blocks = npao_species(the_species)*n_pts_in_block
|
||||
temp_block_storage = zero
|
||||
offset_position = no_of_ib_ia
|
||||
!$omp parallel do default(none) &
|
||||
!$omp schedule(dynamic) &
|
||||
!$omp reduction(+: temp_block_storage) &
|
||||
!$omp shared(n_pts_in_block, direction, pao, the_species, offset_position, &
|
||||
!$omp npoint, no_of_ib_ia, ip_store, r_store, x_store, y_store, z_store) &
|
||||
!$omp private(r_from_i, ip, position, ipoint, &
|
||||
!$omp x, y, z, l1, acz, m1, count1, val)
|
||||
do ip=1,npoint
|
||||
ipoint=ip_store(ip)
|
||||
position= offset_position + ipoint
|
||||
|
|
@ -395,19 +396,23 @@ contains
|
|||
! inside the loop might not be great for vectorization, but with the function
|
||||
! call here, it probably won't vectorize anyways.
|
||||
! This loop nest really should be inside pao_elem_derivative_2, see issue #198
|
||||
gridfunctions(pao_fns)%griddata(position+(count1-1)*n_pts_in_block) = val
|
||||
!gridfunctions(pao_fns)%griddata(position+(count1-1)*n_pts_in_block) = val
|
||||
temp_block_storage(ipoint + (count1-1)*n_pts_in_block) = val
|
||||
count1 = count1+1
|
||||
end do ! m1
|
||||
end do ! acz
|
||||
end do ! l1
|
||||
enddo ! ip=1,npoint
|
||||
!$omp end parallel do
|
||||
gridfunctions(pao_fns)%griddata(no_of_ib_ia+1:no_of_ib_ia+current_num_blocks) = &
|
||||
temp_block_storage(1:current_num_blocks)
|
||||
endif! (npoint > 0) then
|
||||
no_of_ib_ia = no_of_ib_ia + npao_species(the_species)*n_pts_in_block
|
||||
enddo ! naba_atoms
|
||||
enddo ! naba_part
|
||||
endif !(naba_atoms_of_blocks(atomf)%no_of_part(iblock) > 0) !naba atoms?
|
||||
enddo ! iblock : primary set of blocks
|
||||
!$omp end parallel do
|
||||
deallocate(temp_block_storage)
|
||||
call my_barrier()
|
||||
call start_timer(tmr_std_allocation)
|
||||
deallocate(ip_store,x_store,y_store,z_store,r_store)
|
||||
|
|
|
|||
|
|
@ -1207,8 +1207,8 @@ contains
|
|||
|
||||
! 2019Dec30 tsuyoshi
|
||||
! Dump Kmatrix every n_dumpSCF, if n_dumpSCF > 0
|
||||
if(n_dumpSCF > 0 .and. mod(n_iters,n_dumpSCF)==0) then
|
||||
call dump_pos_and_matrices(index=unit_SCF_save)
|
||||
if(n_dumpSCF > 0) then
|
||||
if(mod(n_iters,n_dumpSCF)==0) call dump_pos_and_matrices(index=unit_SCF_save)
|
||||
endif
|
||||
|
||||
! increment iteration counter
|
||||
|
|
|
|||
|
|
@ -229,6 +229,7 @@ contains
|
|||
! Local variables
|
||||
real(double) :: loc_x_energy, exx_tmp
|
||||
|
||||
XC_GGA_stress = zero
|
||||
select case(flag_functional_type)
|
||||
case (functional_lda_pz81)
|
||||
! NOT SPIN POLARISED
|
||||
|
|
|
|||
|
|
@ -366,6 +366,7 @@ contains
|
|||
! Local variables
|
||||
real(double) :: loc_x_energy, exx_tmp
|
||||
|
||||
XC_GGA_stress = zero
|
||||
if(flag_use_libxc) then
|
||||
call get_libxc_potential(density=density, size=size,&
|
||||
xc_potential =xc_potential, &
|
||||
|
|
|
|||
|
|
@ -370,6 +370,7 @@ contains
|
|||
! Local variables
|
||||
real(double) :: loc_x_energy, exx_tmp
|
||||
|
||||
XC_GGA_stress = zero
|
||||
if(flag_use_libxc) then
|
||||
call get_libxc_potential(density=density, size=size,&
|
||||
xc_potential =xc_potential, &
|
||||
|
|
|
|||
|
|
@ -249,7 +249,6 @@ contains
|
|||
send_array(ii+1:ii+nonef)=send_array(ii+1:ii+nonef)+acc_block(i_acc_b+1:i_acc_b+nonef)
|
||||
ii = ii + nonef
|
||||
end do
|
||||
|
||||
Endif ! if the two atoms are within a range
|
||||
enddo ! Loop over right naba atoms
|
||||
enddo ! Loop over left naba atoms
|
||||
|
|
@ -441,7 +440,7 @@ contains
|
|||
use naba_blk_module, only:naba_atm_of_blk
|
||||
use set_blipgrid_module, only: naba_atoms_of_blocks
|
||||
use bucket_module, only:local_bucket,remote_bucket
|
||||
use GenBlas, only:axpy
|
||||
use GenBlas, only:axpy, gemm
|
||||
use comm_array_module, only:send_array
|
||||
use block_module, only:n_pts_in_block
|
||||
use functions_on_grid, only: gridfunctions, fn_on_grid
|
||||
|
|
@ -467,8 +466,6 @@ contains
|
|||
integer :: nsf1, nsf2, ii, stat
|
||||
real(double) :: factor_M
|
||||
|
||||
real(double), dimension(n_pts_in_block) :: griddata_temp_buffer
|
||||
|
||||
call start_timer(tmr_std_integration)
|
||||
!(pointer)
|
||||
! This structure is introduced to keep the strict connection
|
||||
|
|
@ -494,7 +491,7 @@ contains
|
|||
!$omp n_pts_in_block, send_array, gridone, gridtwo) &
|
||||
!$omp private(iprim_blk, n_dim_one, n_dim_two, naba1, naba2, bucket, &
|
||||
!$omp ind_halo1, ind_halo2, nonef, ntwof, nsf1, nsf2, &
|
||||
!$omp ii, factor_M, ind1, ind2, griddata_temp_buffer)
|
||||
!$omp ii, factor_M, ind1, ind2)
|
||||
do iprim_blk=1, domain%groups_on_node
|
||||
|
||||
! In the future we have to prepare n_dim_one & n_dim_two
|
||||
|
|
@ -513,35 +510,20 @@ contains
|
|||
! I HAVE TO BE CAREFUL about WHICH ONE IS LEFT
|
||||
do naba1=1, naba_atm1%no_of_atom(iprim_blk) ! left atom
|
||||
ind_halo1 = naba_atm1%list_atom_by_halo(naba1,iprim_blk)
|
||||
nonef = norb(naba_atm1, naba1, iprim_blk)
|
||||
ind1=n_pts_in_block*(naba_atm1%ibegin_blk_orb(iprim_blk)-1+ &
|
||||
naba_atm1%ibeg_orb_atom(naba1,iprim_blk)-1)+1
|
||||
do naba2=1, naba_atm2%no_of_atom(iprim_blk) ! right atom
|
||||
ind_halo2 = naba_atm2%list_atom_by_halo(naba2,iprim_blk)
|
||||
|
||||
bucket = loc_bucket%i_h2d(ind_halo2,ind_halo1) !index of the pair
|
||||
If(bucket /= 0) then ! naba1 and naba2 makes pair
|
||||
nonef = norb(naba_atm1, naba1, iprim_blk)
|
||||
ntwof = norb(naba_atm2, naba2, iprim_blk)
|
||||
do nsf2=1, ntwof
|
||||
ind2=n_pts_in_block*(naba_atm2%ibegin_blk_orb(iprim_blk)-1+ &
|
||||
naba_atm2%ibeg_orb_atom(naba2,iprim_blk)-1+(nsf2-1))+1
|
||||
do nsf1=1, nonef
|
||||
ind1=n_pts_in_block*(naba_atm1%ibegin_blk_orb(iprim_blk)-1+ &
|
||||
naba_atm1%ibeg_orb_atom(naba1,iprim_blk)-1+(nsf1-1))+1
|
||||
ii = (bucket - 1) + (nsf2 - 1) * nonef + nsf1
|
||||
factor_M=send_array(ii)
|
||||
#ifdef DEBUG
|
||||
call axpy(n_pts_in_block, factor_M, &
|
||||
gridfunctions(gridtwo)%griddata(ind2:ind2+n_pts_in_block-1), 1, &
|
||||
griddata_temp_buffer, 1)
|
||||
!$omp critical
|
||||
gridfunctions(gridone)%griddata(ind1:ind1+n_pts_in_block-1) = griddata_temp_buffer
|
||||
!$omp end critical
|
||||
#else
|
||||
call axpy(n_pts_in_block, factor_M, &
|
||||
gridfunctions(gridtwo)%griddata(ind2:ind2+n_pts_in_block-1), 1, &
|
||||
gridfunctions(gridone)%griddata(ind1:ind1+n_pts_in_block-1), 1)
|
||||
#endif
|
||||
end do
|
||||
end do
|
||||
ind2=n_pts_in_block*(naba_atm2%ibegin_blk_orb(iprim_blk)-1+ &
|
||||
naba_atm2%ibeg_orb_atom(naba2,iprim_blk)-1)+1
|
||||
call gemm('N','T',n_pts_in_block,nonef,ntwof, &
|
||||
one,gridfunctions(gridtwo)%griddata(ind2:ind2+ntwof*n_pts_in_block-1),n_pts_in_block,&
|
||||
send_array(bucket:bucket+nonef*ntwof-1),nonef,&
|
||||
one,gridfunctions(gridone)%griddata(ind1:ind1+nonef*n_pts_in_block-1),n_pts_in_block)
|
||||
Endif ! if the two atoms are within a range
|
||||
|
||||
enddo ! Loop over right naba atoms
|
||||
|
|
|
|||
|
|
@ -106,6 +106,9 @@ module force_module
|
|||
PP_stress, GPV_stress, XC_stress, &
|
||||
nonSCF_stress, pcc_stress, NA_stress
|
||||
|
||||
! How much we mix input and output densities for the calculations of XC_GGA_stress
|
||||
! during non-SCF simulations (defaults to half for averaging)
|
||||
real(double) :: mix_input_output_XC_GGA_stress
|
||||
! Useful parameters for selecting force calculations in NL part
|
||||
integer, parameter :: HF = 1
|
||||
integer, parameter :: Pulay = 2
|
||||
|
|
@ -391,16 +394,15 @@ contains
|
|||
GPV_stress(dir1,dir1) = (hartree_energy_total_rho + &
|
||||
local_ps_energy - core_correction) ! core contains 1/V term
|
||||
end if
|
||||
! XC_GGA_stress is zero for LDA
|
||||
if(flag_self_consistent .or. flag_mix_L_SC_min) then
|
||||
XC_stress(dir1,dir1) = xc_energy + &
|
||||
spin_factor*XC_GGA_stress(dir1,dir1)
|
||||
else ! nonSCF XC found later, along with corrections to Hartree
|
||||
XC_stress(dir1,dir1) = delta_E_xc !xc_energy + spin_factor*XC_GGA_stress(direction)
|
||||
XC_stress(dir1,dir1) = xc_energy + spin_factor*XC_GGA_stress(dir1,dir1)
|
||||
else ! The rest of nonSCF XC found later (nonSC or PCC routines)
|
||||
XC_stress(dir1,dir1) = delta_E_xc
|
||||
end if
|
||||
end do
|
||||
end if
|
||||
WhichPulay = BothPulay
|
||||
|
||||
! matK->matKatomf backtransformation for contracted SFs
|
||||
if (atomf.ne.sf) then
|
||||
do ispin = 1, nspin
|
||||
|
|
@ -761,7 +763,6 @@ contains
|
|||
volume = rcellx*rcelly*rcellz
|
||||
! We need pressure in GPa, and only diagonal terms output
|
||||
scale = -HaBohr3ToGPa/volume
|
||||
!call print_stress("Total pressure: ", stress*scale, 0)
|
||||
if(inode==ionode.AND.iprint_MD + min_layer>=1) &
|
||||
write(io_lun,'(4x,a,f15.8,a4/)') trim(prefix)//" Average pressure: ", &
|
||||
third*scale*(stress(1,1) + stress(2,2) + stress(3,3))," GPa"
|
||||
|
|
@ -3412,9 +3413,13 @@ contains
|
|||
!! Added atomic stress contributions
|
||||
!! 2021/03/23 15:52 dave
|
||||
!! Bugfix: correct calculation of rsq implemented, variables tidied
|
||||
!! 2023/07/19 10:40 dave
|
||||
!! Added minus sign to non-SCF PCC force for consistency with non-SCF part
|
||||
!! 2023/07/24 16:00 dave
|
||||
!! Bug fix for non-SCF NA stress (remove Hartree stress terms)
|
||||
!! SOURCE
|
||||
!!
|
||||
subroutine get_nonSC_correction_force(HF_force, density_out, inode, &
|
||||
subroutine get_nonSC_correction_force(NSCforce, density_out, inode, &
|
||||
ionode, n_atoms, nsize)
|
||||
|
||||
use datatypes
|
||||
|
|
@ -3426,7 +3431,7 @@ contains
|
|||
area_moveatoms, IPRINT_TIME_THRES3, &
|
||||
flag_pcc_global, nspin, spin_factor, &
|
||||
flag_full_stress, flag_stress, &
|
||||
flag_atomic_stress
|
||||
flag_atomic_stress, flag_neutral_atom
|
||||
use XC, only: get_xc_potential, &
|
||||
get_dxc_potential, &
|
||||
flag_is_GGA
|
||||
|
|
@ -3440,7 +3445,7 @@ contains
|
|||
use atomic_density, only: atomic_density_table
|
||||
use pseudo_tm_info, only: pseudo
|
||||
use dimens, only: grid_point_volume, n_my_grid_points
|
||||
use GenBlas, only: axpy
|
||||
use GenBlas, only: axpy, dot
|
||||
use density_module, only: density, density_scale, density_pcc
|
||||
use hartree_module, only: hartree, Hartree_stress
|
||||
use potential_module, only: potential
|
||||
|
|
@ -3451,6 +3456,7 @@ contains
|
|||
WITH_LEVEL, TIME_ACCUMULATE_NO, &
|
||||
TIME_ACCUMULATE_YES
|
||||
use pseudopotential_common, only: pseudopotential
|
||||
use XC, ONLY: XC_GGA_stress
|
||||
|
||||
implicit none
|
||||
|
||||
|
|
@ -3458,7 +3464,7 @@ contains
|
|||
integer :: n_atoms, nsize
|
||||
integer :: inode, ionode
|
||||
real(double), dimension(:,:) :: density_out
|
||||
real(double), dimension(:,:) :: HF_force
|
||||
real(double), dimension(:,:) :: NSCforce
|
||||
|
||||
! Local variables
|
||||
character(len=80) :: sub_name = "get_nonSC_correction_force"
|
||||
|
|
@ -3481,17 +3487,11 @@ contains
|
|||
type(cq_timer) :: backtrace_timer
|
||||
|
||||
real(double), dimension(3,3) :: loc_stress
|
||||
real(double), dimension(:), allocatable :: h_potential, &
|
||||
density_total, &
|
||||
density_out_total
|
||||
real(double), dimension(:), allocatable :: h_potential, density_total, density_out_total
|
||||
real(double), dimension(:,:,:), allocatable :: dVxc_drho
|
||||
real(double), dimension(nspin) :: pot_here, pot_here_pcc
|
||||
! only for GGA with P.C.C.
|
||||
real(double), allocatable, dimension(:) :: h_potential_in, &
|
||||
wk_grid_total, &
|
||||
density_out_GGA_total
|
||||
real(double), allocatable, dimension(:,:) :: wk_grid, &
|
||||
density_out_GGA
|
||||
real(double), allocatable, dimension(:,:) :: wk_grid, density_out_GGA
|
||||
|
||||
|
||||
!****lat<$
|
||||
|
|
@ -3500,36 +3500,24 @@ contains
|
|||
! Spin-polarised PBE non-SCF forces not implemented, so exit if necessary
|
||||
if ((nspin == 2) .and. flag_is_GGA) then ! Only true for CQ not LibXC
|
||||
call cq_warn(sub_name, "NonSCF forces not implemented for spin and GGA; these will be set to zero.")
|
||||
HF_force = zero
|
||||
NSCforce = zero
|
||||
return
|
||||
end if
|
||||
|
||||
stat = 0
|
||||
call start_timer(tmr_std_allocation)
|
||||
allocate(h_potential(nsize), STAT=stat)
|
||||
if (stat /= 0) &
|
||||
call cq_abort("get_nonSC_correction_force: Error alloc mem: ", nsize)
|
||||
if (stat /= 0) call cq_abort("get_nonSC_correction_force: Error alloc mem: ", nsize)
|
||||
allocate(density_total(nsize), STAT=stat)
|
||||
if (stat /= 0) &
|
||||
call cq_abort("get_nonSC_correction_force: Error alloc mem: ", nsize)
|
||||
if (stat /= 0) call cq_abort("get_nonSC_correction_force: Error alloc mem: ", nsize)
|
||||
allocate(density_out_total(nsize), STAT=stat)
|
||||
if (stat /= 0) &
|
||||
call cq_abort("get_nonSC_correction_force: Error alloc mem: ", nsize)
|
||||
if (stat /= 0) call cq_abort("get_nonSC_correction_force: Error alloc mem: ", nsize)
|
||||
allocate(dVxc_drho(nsize,nspin,nspin), STAT=stat)
|
||||
if (stat /= 0) &
|
||||
call cq_abort("get_nonSC_correction_force: Error alloc mem: ", nsize)
|
||||
if (stat /= 0) call cq_abort("get_nonSC_correction_force: Error alloc mem: ", nsize)
|
||||
call reg_alloc_mem(area_moveatoms, (3+nspin*nspin)*nsize, type_dbl)
|
||||
call stop_timer(tmr_std_allocation)
|
||||
|
||||
!****lat<$
|
||||
!print*, size(density_total, dim=1)
|
||||
!print*, size(dVxc_drho, dim=1)
|
||||
!print*, size(density_out_total,dim=1)
|
||||
!print*, size(potential, dim=1)
|
||||
!print*, size(density_out_total,dim=1)
|
||||
!****lat>$
|
||||
|
||||
HF_force = zero
|
||||
NSCforce = zero
|
||||
|
||||
dcellx_block = rcellx / blocks%ngcellx
|
||||
dcelly_block = rcelly / blocks%ngcelly
|
||||
|
|
@ -3561,10 +3549,16 @@ contains
|
|||
jacobian = zero
|
||||
! use density_total (input charge) WITHOUT a factor of half so that this term corrects the GPV
|
||||
! found in the main force routine
|
||||
do ipoint = 1,nsize
|
||||
jacobian = jacobian + (density_out_total(ipoint) - density_total(ipoint))* &
|
||||
(h_potential(ipoint) + pseudopotential(ipoint)) ! Also calculate the correction to GPV for local potential
|
||||
end do
|
||||
jacobian = jacobian + dot(nsize,density_out_total,1,pseudopotential,1) - &
|
||||
dot(nsize,density_total,1,pseudopotential,1)
|
||||
if(flag_neutral_atom) then
|
||||
! These should be zero so ensure that they are; also don't include h_potential in jacobian
|
||||
Hartree_stress = zero
|
||||
loc_stress = zero
|
||||
else
|
||||
jacobian = jacobian + dot(nsize,density_out_total,1,h_potential,1) - &
|
||||
dot(nsize,density_total,1,h_potential,1)
|
||||
end if
|
||||
! Don't apply gsum to jacobian because nonSCF_stress will be summed (end of this routine)
|
||||
jacobian = jacobian*grid_point_volume
|
||||
! loc_stress is \int n^{out} V^{PAD}_{Har}
|
||||
|
|
@ -3584,13 +3578,33 @@ contains
|
|||
! DeltaXC is added in the main force routine
|
||||
! For PCC we will do this in the PCC force routine (easier)
|
||||
if (.NOT.flag_pcc_global) then
|
||||
call get_xc_potential(density, dVxc_drho(:,:,1), &
|
||||
potential(:,1), y_pcc, nsize)
|
||||
! Find XC_GGA_stress for density_out - we will average the input and output
|
||||
! as an approximation to the correct (hideously complex) stress so add factor of half
|
||||
call get_xc_potential(density_out, potential(:,:), &
|
||||
dVxc_drho(:,1,1), h_energy, nsize) ! NB dVxc_drho is a dummy here
|
||||
if (flag_stress) then
|
||||
XC_stress(1,1) = XC_stress(1,1) + spin_factor*XC_GGA_stress(1,1) * &
|
||||
mix_input_output_XC_GGA_stress
|
||||
XC_stress(2,2) = XC_stress(2,2) + spin_factor*XC_GGA_stress(2,2) * &
|
||||
mix_input_output_XC_GGA_stress
|
||||
XC_stress(3,3) = XC_stress(3,3) + spin_factor*XC_GGA_stress(3,3) * &
|
||||
mix_input_output_XC_GGA_stress
|
||||
end if
|
||||
! Find XC potential for input density
|
||||
call get_xc_potential(density, potential(:,:), &
|
||||
dVxc_drho(:,1,1), h_energy, nsize) ! NB dVxc_drho is a dummy here
|
||||
! And now add XC_GGA_stress for density_in scaled by half
|
||||
if (flag_stress) then
|
||||
XC_stress(1,1) = XC_stress(1,1) + spin_factor*XC_GGA_stress(1,1) * &
|
||||
(one - mix_input_output_XC_GGA_stress)
|
||||
XC_stress(2,2) = XC_stress(2,2) + spin_factor*XC_GGA_stress(2,2) * &
|
||||
(one - mix_input_output_XC_GGA_stress)
|
||||
XC_stress(3,3) = XC_stress(3,3) + spin_factor*XC_GGA_stress(3,3) * &
|
||||
(one - mix_input_output_XC_GGA_stress)
|
||||
end if
|
||||
jacobian = zero
|
||||
do spin = 1, nspin
|
||||
do ipoint = 1,nsize
|
||||
jacobian = jacobian + spin_factor*density_out(ipoint,spin)*dVxc_drho(ipoint,spin,1)
|
||||
end do
|
||||
jacobian = jacobian + spin_factor*dot(nsize,density_out(:,spin),1,potential(:,spin),1)
|
||||
end do
|
||||
jacobian = jacobian*grid_point_volume
|
||||
call gsum(jacobian) ! gsum as XC_stress isn't summed elsewhere
|
||||
|
|
@ -3610,46 +3624,34 @@ contains
|
|||
call stop_timer(tmr_l_tmp2, TIME_ACCUMULATE_NO)
|
||||
|
||||
! for P.C.C.
|
||||
call start_timer (tmr_l_tmp1, WITH_LEVEL)
|
||||
if (flag_pcc_global) then
|
||||
allocate(wk_grid_total(nsize), wk_grid(nsize,nspin), STAT=stat)
|
||||
wk_grid_total = zero
|
||||
allocate(wk_grid(nsize,nspin), STAT=stat)
|
||||
wk_grid = zero
|
||||
if (stat /= 0) &
|
||||
call cq_abort('Error allocating wk_grids in &
|
||||
&get_nonSC_correction ', stat)
|
||||
call reg_alloc_mem(area_moveatoms, (nspin + 1) * nsize, type_dbl)
|
||||
call reg_alloc_mem(area_moveatoms, nspin * nsize, type_dbl)
|
||||
do spin = 1, nspin
|
||||
wk_grid(:,spin) = density(:,spin) + half * density_pcc(:)
|
||||
wk_grid_total(:) = wk_grid_total(:) + spin_factor * wk_grid(:,spin)
|
||||
end do
|
||||
! only for GGA
|
||||
! Find dVxc_drho (different for LDA and GGA)
|
||||
if (flag_is_GGA) then
|
||||
allocate(density_out_GGA_total(nsize), density_out_GGA(nsize,nspin), STAT=stat)
|
||||
allocate(density_out_GGA(nsize,nspin), STAT=stat)
|
||||
if (stat /= 0) call cq_abort ('Error allocating &
|
||||
&density_out_GGAs in get_nonSC_force ', stat)
|
||||
call reg_alloc_mem(area_moveatoms, (nspin + 1) * nsize, type_dbl)
|
||||
density_out_GGA_total = zero
|
||||
density_out_GGA = zero
|
||||
call reg_alloc_mem(area_moveatoms, nspin * nsize, type_dbl)
|
||||
density_out_GGA = zero
|
||||
do spin = 1, nspin
|
||||
density_out_GGA(:,spin) = density_out(:,spin) + half * density_pcc(:)
|
||||
density_out_GGA_total(:) = density_out_GGA_total(:) + spin_factor * density_out_GGA(:,spin)
|
||||
end do
|
||||
! copy hartree potential
|
||||
allocate(h_potential_in(nsize), STAT=stat)
|
||||
if (stat /= 0) &
|
||||
call cq_abort('Error allocating h_potential_in in &
|
||||
&get_nonSC_force ', stat)
|
||||
call reg_alloc_mem(area_moveatoms, nsize, type_dbl)
|
||||
h_potential_in = h_potential
|
||||
end if
|
||||
end if
|
||||
|
||||
call start_timer (tmr_l_tmp1, WITH_LEVEL)
|
||||
if (flag_pcc_global) then
|
||||
if(flag_is_GGA) then
|
||||
call get_dxc_potential(wk_grid, dVxc_drho, nsize, density_out_GGA)
|
||||
! GGA with spin not implemented !
|
||||
potential(:,1) = potential(:,1) + dVxc_drho(:,1,1)
|
||||
deallocate(density_out_GGA, STAT=stat)
|
||||
if (stat /= 0) call cq_abort('Error deallocating density_out_GGAs in &
|
||||
&get_nonSC_force ', stat)
|
||||
call reg_dealloc_mem(area_moveatoms, nspin * nsize, type_dbl)
|
||||
else
|
||||
call get_dxc_potential(wk_grid, dVxc_drho, nsize)
|
||||
end if
|
||||
|
|
@ -3662,22 +3664,8 @@ contains
|
|||
call get_dxc_potential(density, dVxc_drho, nsize)
|
||||
end if
|
||||
end if
|
||||
! deallocating density_out_GGA: only for P.C.C.
|
||||
if (flag_pcc_global) then
|
||||
if (flag_is_GGA) then
|
||||
deallocate(density_out_GGA_total, density_out_GGA, STAT=stat)
|
||||
if (stat /= 0) call cq_abort('Error deallocating density_out_GGAs in &
|
||||
&get_nonSC_force ', stat)
|
||||
call reg_dealloc_mem(area_moveatoms, (nspin + 1) * nsize, type_dbl)
|
||||
! make a copy of potential at this point
|
||||
! use wk_grid as a temporary storage
|
||||
do spin = 1, nspin
|
||||
wk_grid(:,spin) = potential(:,spin)
|
||||
end do
|
||||
end if
|
||||
end if !flag_pcc_global
|
||||
|
||||
! for LDA
|
||||
! for LDA: find dn* dxc_potential
|
||||
if (.NOT.(flag_is_GGA)) then
|
||||
do spin = 1, nspin
|
||||
do spin_2 = 1, nspin
|
||||
|
|
@ -3696,12 +3684,12 @@ contains
|
|||
|
||||
! Restart of the timer; level assigned here
|
||||
call start_timer(tmr_l_tmp2, WITH_LEVEL)
|
||||
! Calculate the Hartree potential from the output density
|
||||
h_potential = zero
|
||||
! Preserve the Hartree stress we've calculated
|
||||
loc_stress = Hartree_stress
|
||||
loc_stress = Hartree_stress ! Preserve Hartree stress
|
||||
call hartree(density_out_total, h_potential, nsize, h_energy)
|
||||
! And restore
|
||||
Hartree_stress = loc_stress
|
||||
! And subtract from potential so that we have delta V_Ha
|
||||
do spin = 1, nspin
|
||||
call axpy(nsize, -one, h_potential, 1, potential(:,spin), 1)
|
||||
end do
|
||||
|
|
@ -3807,8 +3795,8 @@ contains
|
|||
! on what we are doing here.
|
||||
do spin = 1, nspin
|
||||
do dir1 = 1, 3
|
||||
HF_force(dir1,ig_atom) = &
|
||||
HF_force(dir1,ig_atom) + spin_factor * &
|
||||
NSCforce(dir1,ig_atom) = &
|
||||
NSCforce(dir1,ig_atom) + spin_factor * &
|
||||
fr_1(dir1,spin) * pot_here(spin)
|
||||
if (flag_stress) then
|
||||
if (flag_full_stress) then
|
||||
|
|
@ -3851,23 +3839,16 @@ contains
|
|||
! for GGA
|
||||
potential = zero
|
||||
do spin = 1, nspin
|
||||
do i = 1, n_my_grid_points
|
||||
! -delta n * dxc_potential
|
||||
potential(i,spin) = wk_grid(i,spin) - h_potential_in(i)
|
||||
end do
|
||||
potential(:,spin) = dVxc_drho(:,spin,spin)
|
||||
end do
|
||||
else
|
||||
! For LDA
|
||||
potential = zero
|
||||
do spin = 1, nspin
|
||||
do spin_2 = 1, nspin
|
||||
do i = 1, n_my_grid_points
|
||||
potential(i,spin) = &
|
||||
potential(i,spin) + &
|
||||
spin_factor * &
|
||||
(density(i,spin_2) - density_out(i,spin_2)) * &
|
||||
dVxc_drho(i,spin_2,spin)
|
||||
end do
|
||||
potential(:,spin) = potential(:,spin) + &
|
||||
spin_factor * (density(:,spin_2) - density_out(:,spin_2)) * &
|
||||
dVxc_drho(:,spin_2,spin)
|
||||
end do
|
||||
end do
|
||||
end if
|
||||
|
|
@ -3962,7 +3943,8 @@ contains
|
|||
! calculated from set_density
|
||||
do spin = 1, nspin
|
||||
do dir1 = 1, 3
|
||||
fr_pcc(dir1,spin) = r_pcc(dir1) * half * derivative_pcc * density_scale(spin)
|
||||
fr_pcc(dir1,spin) = -r_pcc(dir1) * half * &
|
||||
derivative_pcc * density_scale(spin)
|
||||
end do
|
||||
end do
|
||||
end if
|
||||
|
|
@ -3974,8 +3956,8 @@ contains
|
|||
! components)
|
||||
do spin = 1, nspin
|
||||
do dir1 = 1, 3
|
||||
HF_force(dir1,ig_atom) = &
|
||||
HF_force(dir1,ig_atom) + spin_factor * &
|
||||
NSCforce(dir1,ig_atom) = &
|
||||
NSCforce(dir1,ig_atom) + spin_factor * &
|
||||
fr_pcc(dir1,spin) * pot_here_pcc(spin)
|
||||
if (flag_stress) then
|
||||
if (flag_full_stress) then
|
||||
|
|
@ -4011,7 +3993,7 @@ contains
|
|||
end if ! (flag_pcc_global)
|
||||
|
||||
call start_timer(tmr_l_tmp1, WITH_LEVEL)
|
||||
call gsum(HF_force, 3, n_atoms)
|
||||
call gsum(NSCforce, 3, n_atoms)
|
||||
if (flag_stress) call gsum(nonSCF_stress,3,3)
|
||||
call stop_print_timer(tmr_l_tmp1, "NSC force - Compilation", &
|
||||
IPRINT_TIME_THRES3)
|
||||
|
|
@ -4019,18 +4001,11 @@ contains
|
|||
! deallocating temporary arrays
|
||||
call start_timer(tmr_std_allocation)
|
||||
if (flag_pcc_global) then
|
||||
deallocate(wk_grid_total, wk_grid, STAT=stat)
|
||||
deallocate(wk_grid, STAT=stat)
|
||||
if (stat /= 0) &
|
||||
call cq_abort('Error deallocating wk_grid in &
|
||||
&get_nonSC_correction_force ', stat)
|
||||
call reg_dealloc_mem(area_moveatoms, (nspin + 1) * nsize, type_dbl)
|
||||
if (flag_is_GGA) then
|
||||
deallocate(h_potential_in, STAT=stat)
|
||||
if (stat /= 0) &
|
||||
call cq_abort('Error deallocating h_potential_in in &
|
||||
&get_nonSC_correction_force ', stat)
|
||||
call reg_dealloc_mem(area_moveatoms, nsize, type_dbl)
|
||||
end if ! for GGA
|
||||
end if ! flag_pcc_global
|
||||
deallocate(h_potential, density_total, density_out_total, dVxc_drho, &
|
||||
STAT=stat)
|
||||
|
|
@ -4142,6 +4117,7 @@ contains
|
|||
print_timer, stop_print_timer, &
|
||||
WITH_LEVEL, TIME_ACCUMULATE_NO, &
|
||||
TIME_ACCUMULATE_YES
|
||||
use XC, ONLY: XC_GGA_stress
|
||||
|
||||
implicit none
|
||||
|
||||
|
|
@ -4171,15 +4147,14 @@ contains
|
|||
real(double), dimension(nspin) :: pot_here_pcc
|
||||
real(double), dimension(3) :: r_pcc, fr_pcc, r
|
||||
! allocatable arrays
|
||||
real(double), dimension(:), allocatable :: xc_epsilon, density_wk_tot
|
||||
real(double), dimension(:), allocatable :: xc_epsilon
|
||||
real(double), dimension(:,:), allocatable :: xc_potential, density_wk
|
||||
type(cq_timer) :: backtrace_timer
|
||||
|
||||
call start_backtrace(t=backtrace_timer,who='get_PCC_force',where=7,level=3,echo=.true.)
|
||||
allocate(xc_epsilon(size), density_wk_tot(size), &
|
||||
xc_potential(size,nspin), density_wk(size,nspin), STAT=stat)
|
||||
allocate(xc_epsilon(size), xc_potential(size,nspin), density_wk(size,nspin), STAT=stat)
|
||||
if (stat /= 0) call cq_abort("get_pcc_force: Error alloc mem: ", size)
|
||||
call reg_alloc_mem(area_moveatoms, (2+2*nspin)*size, type_dbl)
|
||||
call reg_alloc_mem(area_moveatoms, (1+2*nspin)*size, type_dbl)
|
||||
|
||||
! initialise arrays
|
||||
pcc_force = zero
|
||||
|
|
@ -4195,19 +4170,47 @@ contains
|
|||
dcellz_grid = dcellz_block / nz_in_block
|
||||
|
||||
call start_timer (tmr_l_tmp2)
|
||||
! Find XC_GGA_stress using output density for non-SCF case
|
||||
! we will average the input and output as an approximation
|
||||
! to the correct (hideously complex) stress so add mixing factor
|
||||
if((.not.flag_self_consistent).and.(.not.flag_mix_L_SC_min)) then
|
||||
if(.NOT.present(density_out)) call cq_abort("Output density not passed to PCC force for nonSCF calculation")
|
||||
if (flag_stress) then
|
||||
density_wk = zero
|
||||
do spin = 1, nspin
|
||||
density_wk(:,spin) = density_out(:,spin) + half * density_pcc(:)
|
||||
end do
|
||||
call get_xc_potential(density_wk, xc_potential, &
|
||||
xc_epsilon, xc_energy, size)
|
||||
XC_stress(1,1) = XC_stress(1,1) + spin_factor*XC_GGA_stress(1,1) * &
|
||||
mix_input_output_XC_GGA_stress
|
||||
XC_stress(2,2) = XC_stress(2,2) + spin_factor*XC_GGA_stress(2,2) * &
|
||||
mix_input_output_XC_GGA_stress
|
||||
XC_stress(3,3) = XC_stress(3,3) + spin_factor*XC_GGA_stress(3,3) * &
|
||||
mix_input_output_XC_GGA_stress
|
||||
end if
|
||||
end if
|
||||
density_wk = zero
|
||||
density_wk_tot = zero
|
||||
do spin = 1, nspin
|
||||
density_wk(:,spin) = density(:,spin) + half * density_pcc(:)
|
||||
density_wk_tot(:) = density_wk_tot(:) + spin_factor * density_wk(:,spin)
|
||||
end do
|
||||
|
||||
call get_xc_potential(density_wk, xc_potential, &
|
||||
xc_epsilon, xc_energy, size)
|
||||
! And now add XC_GGA_stress for density_in scaled by half
|
||||
if((.not.flag_self_consistent).and.(.not.flag_mix_L_SC_min)) then
|
||||
if (flag_stress) then
|
||||
XC_stress(1,1) = XC_stress(1,1) + spin_factor*XC_GGA_stress(1,1) * &
|
||||
(one - mix_input_output_XC_GGA_stress)
|
||||
XC_stress(2,2) = XC_stress(2,2) + spin_factor*XC_GGA_stress(2,2) * &
|
||||
(one - mix_input_output_XC_GGA_stress)
|
||||
XC_stress(3,3) = XC_stress(3,3) + spin_factor*XC_GGA_stress(3,3) * &
|
||||
(one - mix_input_output_XC_GGA_stress)
|
||||
end if
|
||||
end if
|
||||
if(PRESENT(xc_energy_ret)) xc_energy_ret = xc_energy
|
||||
! We do this here to re-use xc_potential - for non-PCC we do it in get_nonSC_correction_force
|
||||
if((.not.flag_self_consistent).and.(.not.flag_mix_L_SC_min)) then
|
||||
if(.NOT.present(density_out)) call cq_abort("Output density not passed to PCC force for nonSCF calculation")
|
||||
if (flag_stress) then
|
||||
jacobian = zero
|
||||
do spin=1,nspin
|
||||
|
|
@ -4362,9 +4365,9 @@ contains
|
|||
call stop_print_timer(tmr_l_tmp1, "PCC force - Compilation", &
|
||||
IPRINT_TIME_THRES3)
|
||||
|
||||
deallocate(xc_epsilon, density_wk_tot, xc_potential, density_wk, STAT=stat)
|
||||
deallocate(xc_epsilon, xc_potential, density_wk, STAT=stat)
|
||||
if (stat /= 0) call cq_abort("get_pcc_force: Error dealloc mem")
|
||||
call reg_dealloc_mem(area_moveatoms, (2+2*nspin)*size, type_dbl)
|
||||
call reg_dealloc_mem(area_moveatoms, (1+2*nspin)*size, type_dbl)
|
||||
call stop_backtrace(t=backtrace_timer,who='get_PCC_force',echo=.true.)
|
||||
|
||||
return
|
||||
|
|
|
|||
|
|
@ -941,6 +941,7 @@ contains
|
|||
use biblio, only: flag_dump_bib
|
||||
!2019/12/27 tsuyoshi
|
||||
use density_module, only: method_UpdateChargeDensity,DensityMatrix,AtomicCharge
|
||||
use force_module, only: mix_input_output_XC_GGA_stress
|
||||
|
||||
implicit none
|
||||
|
||||
|
|
@ -1387,7 +1388,9 @@ contains
|
|||
! Added DRB 2018/07/16 for safety
|
||||
if(flag_Multisite) then
|
||||
RadiusMS(i) = fdf_double ('Atom.MultisiteRange',zero)
|
||||
RadiusLD(i) = fdf_double ('Atom.LFDRange',zero)
|
||||
RadiusLD(i) = fdf_double ('Atom.LFDRange',RadiusMS(i))
|
||||
if(RadiusLD(i)<RadiusMS(i)) call cq_warn(sub_name,"LFD range should be larger than MSSF range: ", &
|
||||
RadiusLD(i),RadiusMS(i))
|
||||
end if
|
||||
! Moved to ... so that RadiusAtomf is read from ion files
|
||||
!if (flag_Multisite) RadiusSupport(i) = RadiusAtomf(i) + RadiusMS(i)
|
||||
|
|
@ -1570,6 +1573,10 @@ contains
|
|||
sqnm_trust_step = fdf_double ('AtomMove.MaxSQNMStep',0.2_double )
|
||||
LBFGS_history = fdf_integer('AtomMove.LBFGSHistory', 5 )
|
||||
flag_opt_cell = fdf_boolean('AtomMove.OptCell', .false.)
|
||||
! At present (2023/07/26 just before v1.2 release) neutral atom is required for cell opt
|
||||
if(flag_opt_cell.and.(.not.flag_neutral_atom)) &
|
||||
call cq_abort("You must use neutral atom for cell optimisation")
|
||||
! This can be removed when ewald update is implemented
|
||||
flag_variable_cell = flag_opt_cell
|
||||
optcell_method = fdf_integer('AtomMove.OptCellMethod', 1)
|
||||
cell_constraint_flag = fdf_string(20,'AtomMove.OptCell.Constraint','none')
|
||||
|
|
@ -1581,6 +1588,7 @@ contains
|
|||
flag_stress = fdf_boolean('AtomMove.CalcStress', .true.)
|
||||
flag_full_stress = fdf_boolean('AtomMove.FullStress', .false.)
|
||||
flag_atomic_stress = fdf_boolean('AtomMove.AtomicStress', .false.)
|
||||
mix_input_output_XC_GGA_stress = fdf_double('General.MixXCGGAInOut',half)
|
||||
!
|
||||
flag_vary_basis = fdf_boolean('minE.VaryBasis', .false.)
|
||||
if(.NOT.flag_vary_basis) then
|
||||
|
|
@ -2945,7 +2953,7 @@ contains
|
|||
real(double), dimension(1:3) :: mp_shift
|
||||
real(double), allocatable, dimension(:,:) :: kk_tmp
|
||||
real(double), allocatable, dimension(:) :: wtk_tmp
|
||||
integer :: nkp_tmp
|
||||
integer :: nkp_tmp, nkp_in_line, inc
|
||||
integer :: counter
|
||||
character(len=2) :: suffix
|
||||
|
||||
|
|
@ -3067,41 +3075,69 @@ contains
|
|||
write(io_lun,fmt='(4x,"Using ",i3," lines of k-points specified by user")')
|
||||
end if
|
||||
if(nkp_lines<1) call cq_abort("Need to specify how many kpoint lines !",nkp_lines)
|
||||
nkp = fdf_integer('Diag.NumKpts',2)
|
||||
if(iprint_init>1.AND.inode==ionode) write(io_lun,fmt='(8x,"Number of Kpoints in a line: ",i4)') nkp
|
||||
allocate(kk(3,nkp*nkp_lines),wtk(nkp*nkp_lines),STAT=stat)
|
||||
if(stat/=0) call cq_abort('FindEvals: couldnt allocate kpoints',nkp*nkp_lines)
|
||||
call reg_alloc_mem(area_general,4*nkp*nkp_lines,type_dbl)
|
||||
kk = zero
|
||||
wtk = one/real(nkp*nkp_lines,double)
|
||||
sum = zero
|
||||
nk_st = 1
|
||||
nkp_in_line = fdf_integer('Diag.NumKpts',2)
|
||||
if(iprint_init>1.AND.inode==ionode) write(io_lun,fmt='(8x,"Number of Kpoints in a line: ",i4)') nkp_in_line
|
||||
! Total number of k-points (potentially corrected later)
|
||||
nkp = nkp_in_line*nkp_lines
|
||||
! Read start/end points for lines and correct for duplication
|
||||
allocate(kk_tmp(3,2*nkp_lines))
|
||||
kk_tmp = zero
|
||||
if(fdf_block('Diag.KpointLines'))then
|
||||
if(1+block_end-block_start<2*nkp_lines) &
|
||||
call cq_abort("Kpoint line error: ",1+block_end-block_start,nkp_lines)
|
||||
do i=1,2*nkp_lines,2
|
||||
! Start/end points for line
|
||||
read (unit=input_array(block_start+i-1),fmt=*) &
|
||||
kk(1,nk_st),kk(2,nk_st),kk(3,nk_st)
|
||||
kk_tmp(1,i),kk_tmp(2,i),kk_tmp(3,i)
|
||||
read (unit=input_array(block_start+i),fmt=*) &
|
||||
kk(1,nk_st+nkp-1),kk(2,nk_st+nkp-1),kk(3,nk_st+nkp-1)
|
||||
dkx = (kk(1,nk_st+nkp-1) - kk(1,nk_st))/real(nkp-1,double)
|
||||
dky = (kk(2,nk_st+nkp-1) - kk(2,nk_st))/real(nkp-1,double)
|
||||
dkz = (kk(3,nk_st+nkp-1) - kk(3,nk_st))/real(nkp-1,double)
|
||||
if(iprint_init>1.AND.inode==ionode) &
|
||||
write(io_lun,fmt='(2x,"K-point spacing along line : ",i3,3f7.3)') i,dkx,dky,dkz
|
||||
do j=1,nkp-2
|
||||
kk(1,nk_st+j) = kk(1,nk_st+j-1)+dkx
|
||||
kk(2,nk_st+j) = kk(2,nk_st+j-1)+dky
|
||||
kk(3,nk_st+j) = kk(3,nk_st+j-1)+dkz
|
||||
end do
|
||||
nk_st = nk_st + nkp
|
||||
kk_tmp(1,i+1),kk_tmp(2,i+1),kk_tmp(3,i+1)
|
||||
! If the start of this line duplicates the end of the last, reduce nkp by 1
|
||||
if(i>1) then
|
||||
if(abs(kk_tmp(1,i)-kk_tmp(1,i-1))<RD_ERR .and. &
|
||||
abs(kk_tmp(2,i)-kk_tmp(2,i-1))<RD_ERR .and. &
|
||||
abs(kk_tmp(3,i)-kk_tmp(3,i-1))<RD_ERR) nkp = nkp - 1
|
||||
end if
|
||||
end do
|
||||
else
|
||||
call cq_abort("Must specify a block Diag.KpointLines to have lines of kpoints !")
|
||||
end if
|
||||
nkp = nkp*nkp_lines
|
||||
allocate(kk(3,nkp),wtk(nkp),STAT=stat)
|
||||
if(stat/=0) call cq_abort('FindEvals: couldnt allocate kpoints',nkp)
|
||||
call reg_alloc_mem(area_general,4*nkp,type_dbl)
|
||||
kk = zero
|
||||
nk_st = 1
|
||||
do i=1,2*nkp_lines,2
|
||||
! Spacing for this line in three directions
|
||||
dkx = (kk_tmp(1,i+1) - kk_tmp(1,i))/real(nkp_in_line-1,double)
|
||||
dky = (kk_tmp(2,i+1) - kk_tmp(2,i))/real(nkp_in_line-1,double)
|
||||
dkz = (kk_tmp(3,i+1) - kk_tmp(3,i))/real(nkp_in_line-1,double)
|
||||
if(iprint_init>1.AND.inode==ionode) &
|
||||
write(io_lun,fmt='(2x,"K-point spacing along line : ",i3,3f7.3)') i,dkx,dky,dkz
|
||||
! Number of points in line
|
||||
inc = nkp_in_line
|
||||
if(i<2*nkp_lines-1) then ! If last point is the same as first point of next line, ignore
|
||||
if(abs(kk_tmp(1,i+2)-kk_tmp(1,i+1))<RD_ERR .and. &
|
||||
abs(kk_tmp(2,i+2)-kk_tmp(2,i+1))<RD_ERR .and. &
|
||||
abs(kk_tmp(3,i+2)-kk_tmp(3,i+1))<RD_ERR) then
|
||||
inc = nkp_in_line - 1
|
||||
end if
|
||||
end if
|
||||
! Initial point
|
||||
kk(1,nk_st) = kk_tmp(1,i)
|
||||
kk(2,nk_st) = kk_tmp(2,i)
|
||||
kk(3,nk_st) = kk_tmp(3,i)
|
||||
! Intermediate points
|
||||
do j=1,inc-1
|
||||
kk(1,nk_st+j) = kk(1,nk_st+j-1)+dkx
|
||||
kk(2,nk_st+j) = kk(2,nk_st+j-1)+dky
|
||||
kk(3,nk_st+j) = kk(3,nk_st+j-1)+dkz
|
||||
end do
|
||||
nk_st = nk_st + inc
|
||||
end do
|
||||
deallocate(kk_tmp)
|
||||
wtk = one/real(nkp,double)
|
||||
! Write out fractional k-points
|
||||
if(iprint_init>0.AND.inode==ionode) then
|
||||
if(iprint_init>1.AND.inode==ionode) then
|
||||
write(io_lun,7) nkp
|
||||
do i=1,nkp
|
||||
write(io_lun,fmt='(8x,i5,3f15.6,f12.3)')&
|
||||
|
|
|
|||
|
|
@ -3775,7 +3775,7 @@ contains
|
|||
if (.not.flag_LFD_MD_UseAtomicDensity) call set_atomic_density(.false.)
|
||||
end if
|
||||
! If we have read K and are predicting density from it, then rebuild
|
||||
if(flag_diagonalisation.AND.flag_LmatrixReuse) then
|
||||
if(flag_diagonalisation.AND.flag_LmatrixReuse.AND.flag_self_consistent) then
|
||||
call get_electronic_density(density,electrons,atomfns,H_on_atomfns(1), &
|
||||
inode,ionode,maxngrid)
|
||||
do spin=1,nspin
|
||||
|
|
|
|||
|
|
@ -558,8 +558,8 @@ contains
|
|||
!enddo
|
||||
|
||||
! Write out current SF coefficients every n_dumpSFcoeff, if n_dumpSFcoeff > 0)
|
||||
if (n_dumpSFcoeff > 0 .and. mod(n_iterations,n_dumpSFcoeff) == 0) then
|
||||
call dump_pos_and_matrices(index = unit_MSSF_save)
|
||||
if (n_dumpSFcoeff > 0) then
|
||||
if(mod(n_iterations,n_dumpSFcoeff) == 0) call dump_pos_and_matrices(index = unit_MSSF_save)
|
||||
endif
|
||||
|
||||
! Find change in energy for convergence
|
||||
|
|
@ -1181,8 +1181,8 @@ contains
|
|||
endif
|
||||
|
||||
! Write out current SF coefficients and density matrices with some iprint (in future)
|
||||
if(n_dumpSFcoeff > 0 .and. mod(iter,n_dumpSFcoeff) ==0) then
|
||||
call dump_pos_and_matrices(index=unit_MSSF_save)
|
||||
if(n_dumpSFcoeff > 0 ) then
|
||||
if(mod(iter,n_dumpSFcoeff) ==0) call dump_pos_and_matrices(index=unit_MSSF_save)
|
||||
endif
|
||||
! Go out if converged
|
||||
if (convergence_flag) then
|
||||
|
|
|
|||
|
|
@ -31,9 +31,9 @@ LIBS= $(FFT_LIB) $(XC_LIB) -lscalapack $(BLAS)
|
|||
#XC_COMPFLAGS =
|
||||
|
||||
# LibXC compatibility
|
||||
# Choose LibXC version: v4 or v5
|
||||
XC_LIBRARY = LibXC_v4
|
||||
# XC_LIBRARY = LibXC_v5
|
||||
# Choose LibXC version: v4 (deprecated) or v5/6 (v5 and v6 have the same interface)
|
||||
# XC_LIBRARY = LibXC_v4
|
||||
XC_LIBRARY = LibXC_v5
|
||||
XC_LIB = -lxcf90 -lxc
|
||||
XC_COMPFLAGS = -I/usr/local/include
|
||||
|
||||
|
|
|
|||
|
|
@ -22,11 +22,10 @@ These steps can be run automatically using the script `run_conquest_test.sh` in
|
|||
|
||||
To add new tests
|
||||
|
||||
1. Add input files and a sample output file (run with `IO.Iprint 0` and named `Conquest_out.ref`) in a new subdirectory under [testsuite](./). The naming convention is test directory names start with `test_` followed by a running index with three digits, e.g. `003`.
|
||||
2. Add an entry to the `"test_path"` parameter in [test_check_output.py](./test_check_output.py). The test driver will check for all fields in the output listed in the `"key"` parameter. They are read from the [Conquest_out](test_001_bulk_Si_1proc_Diag/Conquest_out.ref) file by the `read_conquest_out()` function.
|
||||
3. Add it as a new `Run test XXX` step to the [CI workflow](../.github/workflows/makefile.yml)
|
||||
4. *optional* If a new field in the output needs to be checked, these things are required:
|
||||
- Add the logic how to read the line of output to the loop `for line in file.readlines()` in `read_conquest_out()`. The result should be stored in the dictionary `Results`.
|
||||
- Add the new key from the dictionary `Results` to the `"key"` parameter for `test_check_outputs()`.
|
||||
- By default all keys are checked in all tests. If the new key creates errors in the previous tests (e.g. it is not found in their output), add the offending combination(s) to the `xfail_test` logical variable, so that failures in those tests are expected and don't throw errors. See the `"Not-a-real-key"` key as an example.
|
||||
5. *optional* By default the fields are checked to 6 decimal precision. If a custom precision is required, add it to the `custom_precision` dictionary with the corresponding key.
|
||||
1. Add input files and a sample output file (run with `IO.Iprint 0` and named `Conquest_out.ref`) in a new subdirectory under [testsuite](./). The naming convention is test directory names start with `test_` followed by a running index with three digits, e.g. `004`.
|
||||
2. Add a new test to `TestClass` in [`test_check_output.py`](./test_check_output.py). You can use one of the current tests, named `test_XXX` as templates.
|
||||
- Update the directory name passed to `path`
|
||||
- Update the list of parameters in the `@pytest.mark.parametrize` decorator. The `key` parameters correspond to fields in the `Conquest_out` file to be checked against the reference.
|
||||
- If necessary, update the `read_conquest_out()` function to parse a new field from the output.
|
||||
- If necessary, update the `precision()` function to return a custom precision for a given value of the `key` parameter.
|
||||
4. Add it as a new `Run test XXX` step to the [CI workflow](../.github/workflows/makefile.yml)
|
||||
|
|
|
|||
|
|
@ -37,47 +37,55 @@ def results(path, key):
|
|||
|
||||
return (ref_result[key], test_result[key])
|
||||
|
||||
def precision(key='_'):
|
||||
'''
|
||||
Return the relative tolerance used by tests. By default returns 1e-4, but
|
||||
takes a key as an argument and you can match it to return a different precision
|
||||
|
||||
For example:
|
||||
|
||||
if(key == 'Special case'):
|
||||
return 999.9
|
||||
else:
|
||||
return 1e-4
|
||||
'''
|
||||
|
||||
return 1e-4
|
||||
|
||||
@pytest.fixture
|
||||
def directory():
|
||||
def testsuite_directory():
|
||||
'''
|
||||
Return path to testsuite
|
||||
'''
|
||||
return pathlib.Path(__file__).parent.resolve()
|
||||
|
||||
@pytest.fixture
|
||||
def precision():
|
||||
'''
|
||||
Return the relative tolerance used by tests
|
||||
'''
|
||||
class TestClass:
|
||||
@pytest.mark.parametrize("key",['Harris-Foulkes energy',
|
||||
'Max force',
|
||||
'Force residual',
|
||||
'Total stress'])
|
||||
def test_001(self, key, testsuite_directory):
|
||||
|
||||
return 1e-4
|
||||
path = os.path.join(testsuite_directory, "test_001_bulk_Si_1proc_Diag")
|
||||
res = results(path, key)
|
||||
np.testing.assert_allclose(res[0], res[1], rtol = precision(key), verbose = True)
|
||||
|
||||
@pytest.mark.parametrize("key",['Harris-Foulkes energy',
|
||||
'Max force',
|
||||
'Force residual',
|
||||
'Total stress'])
|
||||
def test_001(key, precision, directory):
|
||||
@pytest.mark.parametrize("key", ['Harris-Foulkes energy',
|
||||
'Max force',
|
||||
'Force residual',
|
||||
'Total stress'])
|
||||
def test_002(self, key, testsuite_directory):
|
||||
|
||||
path = os.path.join(directory, "test_001_bulk_Si_1proc_Diag")
|
||||
res = results(path, key)
|
||||
np.testing.assert_allclose(res[0], res[1], rtol = precision, verbose = True)
|
||||
path = os.path.join(testsuite_directory, "test_002_bulk_Si_1proc_OrderN")
|
||||
res = results(path, key)
|
||||
np.testing.assert_allclose(res[0], res[1], rtol = precision(key), verbose = True)
|
||||
|
||||
@pytest.mark.parametrize("key", ['Harris-Foulkes energy',
|
||||
'Max force',
|
||||
'Force residual',
|
||||
'Total stress'])
|
||||
def test_002(key, precision, directory):
|
||||
@pytest.mark.parametrize("key", ['Harris-Foulkes energy',
|
||||
'Max force',
|
||||
'Force residual',
|
||||
'Total polarisation'])
|
||||
def test_003(self, key, testsuite_directory):
|
||||
|
||||
path = os.path.join(directory, "test_002_bulk_Si_1proc_OrderN")
|
||||
res = results(path, key)
|
||||
np.testing.assert_allclose(res[0], res[1], rtol = precision, verbose = True)
|
||||
|
||||
@pytest.mark.parametrize("key", ['Harris-Foulkes energy',
|
||||
'Max force',
|
||||
'Force residual',
|
||||
'Total polarisation'])
|
||||
def test_003(key, precision, directory):
|
||||
|
||||
path = os.path.join(directory, "test_003_bulk_BTO_polarisation")
|
||||
res = results(path, key)
|
||||
np.testing.assert_allclose(res[0], res[1], rtol = precision, verbose = True)
|
||||
path = os.path.join(testsuite_directory, "test_003_bulk_BTO_polarisation")
|
||||
res = results(path, key)
|
||||
np.testing.assert_allclose(res[0], res[1], rtol = precision(key), verbose = True)
|
||||
|
|
|
|||
|
|
@ -9,6 +9,7 @@ LINKFLAGS= -L/usr/local/lib
|
|||
ARFLAGS=
|
||||
|
||||
# Compilation flags
|
||||
# NB for gcc10 you need to add -fallow-argument-mismatch
|
||||
COMPFLAGS= -O3 $(XC_COMPFLAGS)
|
||||
COMPFLAGS_F77= $(COMPFLAGS)
|
||||
|
||||
|
|
@ -30,10 +31,9 @@ LIBS= $(FFT_LIB) $(XC_LIB) -lscalapack $(BLAS)
|
|||
#XC_COMPFLAGS =
|
||||
|
||||
# LibXC compatibility
|
||||
# Choose LibXC version: v2, v3 or v4
|
||||
#XC_LIBRARY = LibXC_v2
|
||||
#XC_LIBRARY = LibXC_v3
|
||||
XC_LIBRARY = LibXC_v4
|
||||
# Choose LibXC version: v4 (deprecated) or v5/6 (v5 and v6 have the same interface)
|
||||
# XC_LIBRARY = LibXC_v4
|
||||
XC_LIBRARY = LibXC_v5
|
||||
XC_LIB = -lxcf90 -lxc
|
||||
XC_COMPFLAGS = -I/usr/local/include
|
||||
|
||||
|
|
|
|||
|
|
@ -63,6 +63,9 @@ program ConvertCharge
|
|||
write(*,fmt='(2x,"Creating projected density of states (pDOS); also generating DOS"/)')
|
||||
call process_dos
|
||||
call process_pdos
|
||||
else if(i_job==8) then ! DOS output
|
||||
write(*,fmt='(2x,"Creating band structure"/)')
|
||||
call process_band_structure
|
||||
end if
|
||||
! Read eigenvalues and calculate weight for bands by kpt
|
||||
!if(flag_only_charge) then
|
||||
|
|
|
|||
|
|
@ -34,8 +34,10 @@ module local
|
|||
integer, parameter :: dx = 1
|
||||
integer, parameter :: cube = 2
|
||||
|
||||
integer :: flag_proc_band_str
|
||||
logical :: flag_only_charge, flag_by_kpoint, flag_wf_range, flag_proc_range, flag_procwf_range_Ef
|
||||
logical :: flag_total_iDOS, flag_write_forces, flag_write_spin_moments, flag_l_resolved, flag_lm_resolved
|
||||
logical :: flag_outputWF_real
|
||||
character(len=80) :: charge_stub
|
||||
|
||||
integer :: i_job ! Job type
|
||||
|
|
@ -52,5 +54,6 @@ module local
|
|||
real(double) :: kT
|
||||
! Flags controlling Methfessel-Paxton approximation to step-function
|
||||
integer :: flag_smear_type, iMethfessel_Paxton
|
||||
|
||||
integer :: n_atoms_pDOS
|
||||
integer, dimension(:), allocatable :: pDOS_atom_index
|
||||
end module local
|
||||
|
|
|
|||
|
|
@ -96,7 +96,7 @@ contains
|
|||
use local, ONLY: nkp, efermi, current, nptsx, nptsy, nptsz, eigenvalues, flag_by_kpoint, &
|
||||
n_bands_total, band_active_kp, flag_proc_range, band_full_to_active, evec_coeff,&
|
||||
E_procwf_min, E_procwf_max, flag_procwf_range_Ef, band_proc_no, n_bands_process, &
|
||||
grid_x, grid_y, grid_z, wtk
|
||||
grid_x, grid_y, grid_z, wtk, flag_outputWF_real
|
||||
use output, ONLY: write_cube
|
||||
use global_module, only : nspin
|
||||
use read, ONLY: read_eigenvalues, read_psi_coeffs
|
||||
|
|
@ -104,6 +104,7 @@ contains
|
|||
use global_module, ONLY: ni_in_cell, atom_coord, species_glob
|
||||
use species_module, ONLY: nsf_species, n_species
|
||||
use angular_coeff_routines, ONLY: set_prefac_real, set_fact, set_prefac
|
||||
use GenComms, ONLY: cq_abort
|
||||
|
||||
implicit none
|
||||
|
||||
|
|
@ -114,6 +115,8 @@ contains
|
|||
complex(double_cplx), dimension(:,:,:), allocatable :: psi
|
||||
integer :: i_atom, i_spec, i_l, i_zeta,i_m,j_atom,j_spec,j_l,j_zeta,j_m
|
||||
integer :: i_band, i_pao, j_pao
|
||||
real(double), parameter :: band_integral_tol = 1e-3_double
|
||||
real(double) :: max_band_integral_deviation, integral_deviation
|
||||
|
||||
! Create arrays needed by Conquest PAO routines
|
||||
call set_fact(8)
|
||||
|
|
@ -126,13 +129,25 @@ contains
|
|||
call read_psi_coeffs("Process")
|
||||
allocate(current(nptsx,nptsy,nptsz))
|
||||
allocate(psi(nptsx,nptsy,nptsz))
|
||||
max_band_integral_deviation = zero
|
||||
|
||||
if (flag_outputWF_real .and. (nkp.ne.1)) &
|
||||
call cq_abort("OutputWF_real is available only for Gamma-point calculations.")
|
||||
|
||||
if(flag_proc_range) then
|
||||
Emin = E_procwf_min
|
||||
Emax = E_procwf_max
|
||||
if(flag_by_kpoint) then
|
||||
write(*,fmt='(4x,"Writing bands at each k-point")')
|
||||
else
|
||||
write(*,fmt='(4x,"Summing over k-points")')
|
||||
end if
|
||||
if(flag_procwf_range_Ef) then
|
||||
Emin = efermi + Emin
|
||||
Emax = efermi + Emax
|
||||
end if
|
||||
write(*,fmt='(4x,"Writing bands between ",e12.4," and ",e12.4,"Ha as specified in input file")') &
|
||||
Emin(1),Emax(1)
|
||||
do ispin=1,nspin
|
||||
if(flag_by_kpoint) then ! Separate bands by k-point
|
||||
do band=1,n_bands_total
|
||||
|
|
@ -144,12 +159,18 @@ contains
|
|||
psi = zero
|
||||
current = zero
|
||||
call pao_to_grid(band_full_to_active(band), kp, ispin, psi)
|
||||
current = psi*conjg(psi)
|
||||
if (.not.flag_outputWF_real) then
|
||||
current = psi*conjg(psi) ! band density
|
||||
else
|
||||
current = real(psi) ! only real-part of WF
|
||||
endif
|
||||
write(ci,'("Band",I0.6,"den_kp",I0.3,"S",I0.1)') band, kp, ispin
|
||||
call write_cube(current,ci)
|
||||
integral = gpv*sum(current)
|
||||
integral_deviation = abs(integral - one)
|
||||
max_band_integral_deviation = max(integral_deviation, max_band_integral_deviation)
|
||||
! Check for problems with band integral
|
||||
if(abs(integral - one)>1e-4_double) &
|
||||
if(integral_deviation>band_integral_tol) &
|
||||
write(*,fmt='(4x,"Integral of band ",i5," with energy ",f17.10," is ",f17.10)') &
|
||||
band,eigenvalues(band,kp,ispin),integral
|
||||
end if
|
||||
|
|
@ -165,27 +186,41 @@ contains
|
|||
band_active_kp(band,kp,ispin)==1) then
|
||||
call pao_to_grid(band_full_to_active(band), kp, ispin, psi)
|
||||
integral = gpv*sum(psi*conjg(psi))
|
||||
integral_deviation = abs(integral - one)
|
||||
max_band_integral_deviation = max(integral_deviation, max_band_integral_deviation)
|
||||
! Check for problems with band integral
|
||||
if(abs(integral - one)>1e-4_double) &
|
||||
if(integral_deviation>band_integral_tol) &
|
||||
write(*,fmt='(4x,"Integral of band ",i5," at kp ",i5," is ",f17.10)') &
|
||||
band,kp,integral
|
||||
current = current + psi*conjg(psi)*wtk(kp)
|
||||
if (.not.flag_outputWF_real) then
|
||||
current = current + psi*conjg(psi)*wtk(kp) ! band density
|
||||
else
|
||||
current = current + real(psi)*wtk(kp) ! only real-part of WF
|
||||
endif
|
||||
idum1 = 1
|
||||
end if
|
||||
end do ! kp
|
||||
if(idum1==1) then
|
||||
write(ci,'("Band",I0.6,"den_totS",I0.1)') band, ispin
|
||||
call write_cube(current,ci)
|
||||
integral = gpv*sum(current)
|
||||
integral_deviation = abs(integral - one)
|
||||
max_band_integral_deviation = max(integral_deviation, max_band_integral_deviation)
|
||||
! Check for problems with band integral
|
||||
if(integral_deviation>band_integral_tol) &
|
||||
write(*,fmt='(4x,"Integral of band ",i5," is ",f17.10)') &
|
||||
band,integral
|
||||
end if
|
||||
integral = gpv*sum(current)
|
||||
! Check for problems with band integral
|
||||
if(abs(integral - one)>1e-4_double) &
|
||||
write(*,fmt='(4x,"Integral of band ",i5," is ",f17.10)') &
|
||||
band,integral
|
||||
end do ! bands
|
||||
end if
|
||||
end do
|
||||
else ! User has provided list of bands
|
||||
write(*,fmt='(4x,"Writing ",i4," bands specified in input file")') n_bands_process
|
||||
if(flag_by_kpoint) then
|
||||
write(*,fmt='(4x,"Writing bands at each k-point")')
|
||||
else
|
||||
write(*,fmt='(4x,"Summing over k-points")')
|
||||
end if
|
||||
do ispin=1,nspin
|
||||
if(flag_by_kpoint) then ! Separate bands by k-point
|
||||
do band=1,n_bands_process
|
||||
|
|
@ -196,12 +231,18 @@ contains
|
|||
psi = zero
|
||||
current = zero
|
||||
call pao_to_grid(band_full_to_active(band_proc_no(band)), kp, ispin, psi)
|
||||
current = psi*conjg(psi)
|
||||
if (.not.flag_outputWF_real) then
|
||||
current = psi*conjg(psi) ! band density
|
||||
else
|
||||
current = real(psi) ! only real-part of WF
|
||||
endif
|
||||
write(ci,'("Band",I0.6,"den_kp",I0.3,"S",I0.1)') band_proc_no(band), kp, ispin
|
||||
call write_cube(current,ci)
|
||||
integral = gpv*sum(current)
|
||||
integral_deviation = abs(integral - one)
|
||||
max_band_integral_deviation = max(integral_deviation, max_band_integral_deviation)
|
||||
! Check for problems with band integral
|
||||
if(abs(integral - one)>1e-4_double) &
|
||||
if(integral_deviation>band_integral_tol) &
|
||||
write(*,fmt='(4x,"Integral of psi squared ",i5," with energy ",f17.10," is ",f17.10)') &
|
||||
band,eigenvalues(band,kp,ispin),integral
|
||||
end if
|
||||
|
|
@ -215,24 +256,33 @@ contains
|
|||
if(band_active_kp(band_proc_no(band),kp,ispin)==1) then
|
||||
call pao_to_grid(band_full_to_active(band_proc_no(band)), kp, ispin, psi)
|
||||
integral = gpv*sum(psi*conjg(psi))
|
||||
integral_deviation = abs(integral - one)
|
||||
max_band_integral_deviation = max(integral_deviation, max_band_integral_deviation)
|
||||
! Check for problems with band integral
|
||||
if(abs(integral - one)>1e-4_double) &
|
||||
if(integral_deviation>band_integral_tol) &
|
||||
write(*,fmt='(4x,"Integral of band ",i5," at kp ",i5," is ",f17.10)') &
|
||||
band,kp,integral
|
||||
current = current + psi*conjg(psi)*wtk(kp)
|
||||
if (.not.flag_outputWF_real) then
|
||||
current = current + psi*conjg(psi)*wtk(kp) ! band density
|
||||
else
|
||||
current = current + real(psi)*wtk(kp) ! only real-part of WF
|
||||
endif
|
||||
end if
|
||||
end do ! kp
|
||||
write(ci,'("Band",I0.6,"den_totS",I0.1)') band_proc_no(band), ispin
|
||||
call write_cube(current,ci)
|
||||
integral = gpv*sum(current)
|
||||
integral_deviation = abs(integral - one)
|
||||
max_band_integral_deviation = max(integral_deviation, max_band_integral_deviation)
|
||||
! Check for problems with band integral
|
||||
if(abs(integral - one)>1e-4_double) &
|
||||
if(integral_deviation>band_integral_tol) &
|
||||
write(*,fmt='(4x,"Integral of band ",i5," is ",f17.10)') &
|
||||
band,integral
|
||||
end do ! bands
|
||||
end if
|
||||
end do
|
||||
end if
|
||||
write(*,fmt='(4x,"Largest deviation of band integral from one is ",f8.5)') max_band_integral_deviation
|
||||
return
|
||||
end subroutine process_bands
|
||||
|
||||
|
|
@ -396,7 +446,8 @@ contains
|
|||
use datatypes
|
||||
use numbers, ONLY: zero, RD_ERR, twopi, half, one, two, four, six
|
||||
use local, ONLY: eigenvalues, n_bands_total, nkp, wtk, efermi, flag_total_iDOS, &
|
||||
evec_coeff, scaled_evec_coeff, flag_procwf_range_Ef, flag_l_resolved, flag_lm_resolved, band_full_to_active
|
||||
evec_coeff, scaled_evec_coeff, flag_procwf_range_Ef, flag_l_resolved, flag_lm_resolved, &
|
||||
band_full_to_active, n_atoms_pDOS, pDOS_atom_index
|
||||
use read, ONLY: read_eigenvalues, read_psi_coeffs, read_nprocs_from_blocks
|
||||
use global_module, ONLY: nspin, n_DOS, E_DOS_min, E_DOS_max, sigma_DOS, ni_in_cell, species_glob
|
||||
use units, ONLY: HaToeV
|
||||
|
|
@ -424,6 +475,15 @@ contains
|
|||
else if(flag_l_resolved) then
|
||||
write(*,fmt='(4x,"Resolving in l")')
|
||||
end if
|
||||
if(n_atoms_pDOS==ni_in_cell) then
|
||||
write(*,fmt='(4x,"Producing pDOS for all atoms in cell")')
|
||||
else
|
||||
if(n_atoms_pDOS==1) then
|
||||
write(*,fmt='(4x,"Producing pDOS for ",i7," atom in cell")') n_atoms_pDOS
|
||||
else
|
||||
write(*,fmt='(4x,"Producing pDOS for ",i7," atoms in cell")') n_atoms_pDOS
|
||||
end if
|
||||
end if
|
||||
call read_nprocs_from_blocks
|
||||
if(nspin==1) then
|
||||
spin_fac = two
|
||||
|
|
@ -444,25 +504,25 @@ contains
|
|||
allocate(occ(n_bands_total,nkp))
|
||||
! Set up storage based on pDOS per atom, or l/lm resolved per atom
|
||||
if(flag_lm_resolved) then
|
||||
allocate(pDOS_lm(-max_l:max_l,0:max_l,ni_in_cell,n_DOS,nspin))
|
||||
allocate(pDOS_lm(-max_l:max_l,0:max_l,n_atoms_pDOS,n_DOS,nspin))
|
||||
pDOS_lm = zero
|
||||
allocate(total_electrons_l(0:max_l,ni_in_cell,nspin))
|
||||
allocate(total_electrons_l(0:max_l,n_atoms_pDOS,nspin))
|
||||
total_electrons_l = zero
|
||||
! For total pDOS
|
||||
allocate(pDOS(ni_in_cell,n_DOS,nspin))
|
||||
allocate(pDOS(n_atoms_pDOS,n_DOS,nspin))
|
||||
pDOS = zero
|
||||
else if(flag_l_resolved) then
|
||||
allocate(pDOS_l(0:max_l,ni_in_cell,n_DOS,nspin))
|
||||
allocate(pDOS_l(0:max_l,n_atoms_pDOS,n_DOS,nspin))
|
||||
pDOS_l = zero
|
||||
allocate(total_electrons_l(0:max_l,ni_in_cell,nspin))
|
||||
allocate(total_electrons_l(0:max_l,n_atoms_pDOS,nspin))
|
||||
total_electrons_l = zero
|
||||
! For total pDOS
|
||||
allocate(pDOS(ni_in_cell,n_DOS,nspin))
|
||||
allocate(pDOS(n_atoms_pDOS,n_DOS,nspin))
|
||||
pDOS = zero
|
||||
else
|
||||
allocate(pDOS(ni_in_cell,n_DOS,nspin))
|
||||
allocate(pDOS(n_atoms_pDOS,n_DOS,nspin))
|
||||
pDOS = zero
|
||||
allocate(total_electrons(ni_in_cell,nspin))
|
||||
allocate(total_electrons(n_atoms_pDOS,nspin))
|
||||
total_electrons = zero
|
||||
end if
|
||||
! E_DOS_min and max and sigma_DOS already set in process_dos
|
||||
|
|
@ -491,20 +551,25 @@ contains
|
|||
do i = n_min, n_max
|
||||
Ebin = real(i-1,double)*dE_DOS + E_DOS_min
|
||||
a = (Ebin-eigenvalues(i_band, i_kp, i_spin))/sigma_DOS
|
||||
do i_atom = 1, ni_in_cell
|
||||
i_spec = species_glob(i_atom)
|
||||
do i_atom = 1, n_atoms_pDOS
|
||||
i_spec = species_glob(pDOS_atom_index(i_atom))
|
||||
if(flag_l_resolved .and. flag_lm_resolved) then
|
||||
sf_offset = 0
|
||||
do i_l = 0, pao(i_spec)%greatest_angmom
|
||||
nzeta = pao(i_spec)%angmom(i_l)%n_zeta_in_angmom
|
||||
norbs = nzeta
|
||||
do i_m = -i_l,i_l
|
||||
coeff = zdotc(norbs, evec_coeff(sf_offset+1:sf_offset+norbs,i_atom,i_band_c,i_kp,i_spin),1, &
|
||||
scaled_evec_coeff(sf_offset+1:sf_offset+norbs,i_atom,i_band_c,i_kp,i_spin),1)
|
||||
pDOS_lm(i_m,i_l,i_atom,i,i_spin) = pDOS_lm(i_m,i_l,i_atom,i,i_spin) + &
|
||||
coeff = zdotc(norbs, evec_coeff(sf_offset+1:sf_offset+norbs,pDOS_atom_index(i_atom), &
|
||||
i_band_c,i_kp,i_spin),1, &
|
||||
scaled_evec_coeff(sf_offset+1:sf_offset+norbs,pDOS_atom_index(i_atom), &
|
||||
i_band_c,i_kp,i_spin),1)
|
||||
pDOS_lm(i_m,i_l,i_atom,i,i_spin) = &
|
||||
pDOS_lm(i_m,i_l,i_atom,i,i_spin) + &
|
||||
wtk(i_kp)*pf_DOS*exp(-half*a*a)*coeff
|
||||
pDOS(i_atom,i,i_spin) = pDOS(i_atom,i,i_spin) + wtk(i_kp)*pf_DOS*exp(-half*a*a)*coeff
|
||||
total_electrons_l(i_l,i_atom, i_spin) = total_electrons_l(i_l,i_atom, i_spin) + &
|
||||
pDOS(i_atom,i,i_spin) = pDOS(i_atom,i,i_spin) + &
|
||||
wtk(i_kp)*pf_DOS*exp(-half*a*a)*coeff
|
||||
total_electrons_l(i_l,i_atom, i_spin) = &
|
||||
total_electrons_l(i_l,i_atom, i_spin) + &
|
||||
occ(i_band,i_kp)*wtk(i_kp)*pf_DOS*exp(-half*a*a)*coeff
|
||||
sf_offset = sf_offset + norbs
|
||||
end do
|
||||
|
|
@ -514,17 +579,22 @@ contains
|
|||
do i_l = 0, pao(i_spec)%greatest_angmom
|
||||
nzeta = pao(i_spec)%angmom(i_l)%n_zeta_in_angmom
|
||||
norbs = nzeta*(2*i_l+1)
|
||||
coeff = zdotc(norbs, evec_coeff(sf_offset+1:sf_offset+norbs,i_atom,i_band_c,i_kp,i_spin),1, &
|
||||
scaled_evec_coeff(sf_offset+1:sf_offset+norbs,i_atom,i_band_c,i_kp,i_spin),1)
|
||||
pDOS_l(i_l,i_atom,i,i_spin) = pDOS_l(i_l,i_atom,i,i_spin) + wtk(i_kp)*pf_DOS*exp(-half*a*a)*coeff
|
||||
coeff = zdotc(norbs, evec_coeff(sf_offset+1:sf_offset+norbs,pDOS_atom_index(i_atom), &
|
||||
i_band_c,i_kp,i_spin),1, &
|
||||
scaled_evec_coeff(sf_offset+1:sf_offset+norbs,pDOS_atom_index(i_atom), &
|
||||
i_band_c,i_kp,i_spin),1)
|
||||
pDOS_l(i_l,i_atom,i,i_spin) = pDOS_l(i_l,i_atom,i,i_spin) + &
|
||||
wtk(i_kp)*pf_DOS*exp(-half*a*a)*coeff
|
||||
pDOS(i_atom,i,i_spin) = pDOS(i_atom,i,i_spin) + wtk(i_kp)*pf_DOS*exp(-half*a*a)*coeff
|
||||
total_electrons_l(i_l,i_atom, i_spin) = total_electrons_l(i_l,i_atom, i_spin) + &
|
||||
occ(i_band,i_kp)*wtk(i_kp)*pf_DOS*exp(-half*a*a)*coeff
|
||||
sf_offset = sf_offset + norbs
|
||||
end do
|
||||
else
|
||||
coeff = zdotc(npao_species(i_spec),evec_coeff(1:npao_species(i_spec),i_atom,i_band_c,i_kp,i_spin),1, &
|
||||
scaled_evec_coeff(1:npao_species(i_spec),i_atom,i_band_c,i_kp,i_spin),1)
|
||||
coeff = zdotc(npao_species(i_spec),evec_coeff(1:npao_species(i_spec), &
|
||||
pDOS_atom_index(i_atom),i_band_c,i_kp,i_spin),1, &
|
||||
scaled_evec_coeff(1:npao_species(i_spec), &
|
||||
pDOS_atom_index(i_atom),i_band_c,i_kp,i_spin),1)
|
||||
pDOS(i_atom,i,i_spin) = pDOS(i_atom,i,i_spin) + wtk(i_kp)*pf_DOS*exp(-half*a*a)*coeff
|
||||
total_electrons(i_atom, i_spin) = total_electrons(i_atom, i_spin) + &
|
||||
occ(i_band,i_kp)*wtk(i_kp)*pf_DOS*exp(-half*a*a)*coeff
|
||||
|
|
@ -559,8 +629,8 @@ contains
|
|||
write(*,fmt='(4x," Atom Total l=0 l=1 l=2")')
|
||||
write(fmt_DOS,*) max_l + 2 ! Number of columns
|
||||
fmt_DOS = '(4x,i7,x,'//trim(adjustl(fmt_DOS))//'f11.3)'
|
||||
do i_atom = 1, ni_in_cell
|
||||
write(*,fmt=fmt_DOS) i_atom, sum(total_electrons_l(:,i_atom,1)),total_electrons_l(:,i_atom,1)
|
||||
do i_atom = 1, n_atoms_pDOS
|
||||
write(*,fmt=fmt_DOS) pDOS_atom_index(i_atom), sum(total_electrons_l(:,i_atom,1)),total_electrons_l(:,i_atom,1)
|
||||
end do
|
||||
write(*,fmt='(2x,"Integrated pDOS: ",f11.3," electrons")') sum(total_electrons_l)
|
||||
if(flag_lm_resolved) then
|
||||
|
|
@ -577,8 +647,8 @@ contains
|
|||
end if
|
||||
else
|
||||
write(*,fmt='(4x," Atom Electrons")')
|
||||
do i_atom = 1, ni_in_cell
|
||||
write(*,fmt='(4x,i7,x,f11.3)') i_atom, total_electrons(i_atom,1)
|
||||
do i_atom = 1, n_atoms_pDOS
|
||||
write(*,fmt='(4x,i7,x,f11.3)') pDOS_atom_index(i_atom), total_electrons(i_atom,1)
|
||||
end do
|
||||
write(*,fmt='(2x,"Integrated pDOS: ",f11.3," electrons")') sum(total_electrons)
|
||||
end if
|
||||
|
|
@ -589,8 +659,8 @@ contains
|
|||
write(*,fmt='(4x," Atom Spin Up l=0 l=1 l=2 Spin Down l=0 l=1 l=2")')
|
||||
write(fmt_DOS,*) 2*(max_l + 2) ! Number of columns
|
||||
fmt_DOS = '(4x,i7,x,'//trim(adjustl(fmt_DOS))//'f11.3)'
|
||||
do i_atom = 1, ni_in_cell
|
||||
write(*,fmt=fmt_DOS) i_atom, sum(total_electrons_l(:,i_atom,1)),total_electrons_l(:,i_atom,1), &
|
||||
do i_atom = 1, n_atoms_pDOS
|
||||
write(*,fmt=fmt_DOS) pDOS_atom_index(i_atom), sum(total_electrons_l(:,i_atom,1)),total_electrons_l(:,i_atom,1), &
|
||||
sum(total_electrons_l(:,i_atom,2)),total_electrons_l(:,i_atom,2)
|
||||
end do
|
||||
write(*,fmt='(2x,"Integrated spin up pDOS: ",f11.3," electrons")') sum(total_electrons_l(:,:,1))
|
||||
|
|
@ -611,22 +681,22 @@ contains
|
|||
write(*,fmt='(2x,"Results of integrating pDOS between ",f11.3," and ",f11.3," Ha (electrons per atom).")') &
|
||||
E_DOS_min, E_DOS_max
|
||||
write(*,fmt='(4x," Atom Spin Up Spin Down")')
|
||||
do i_atom = 1, ni_in_cell
|
||||
write(*,fmt='(4x,i7,x,2f11.3)') i_atom, total_electrons(i_atom,1), total_electrons(i_atom,2)
|
||||
do i_atom = 1, n_atoms_pDOS
|
||||
write(*,fmt='(4x,i7,x,2f11.3)') pDOS_atom_index(i_atom), total_electrons(i_atom,1), total_electrons(i_atom,2)
|
||||
end do
|
||||
write(*,fmt='(2x,"Integrated spin up pDOS: ",f11.3," electrons")') sum(total_electrons(:,1))
|
||||
write(*,fmt='(2x,"Integrated spin dn pDOS: ",f11.3," electrons")') sum(total_electrons(:,2))
|
||||
end if
|
||||
end if
|
||||
! Write out DOS, shifted to Ef = 0
|
||||
do i_atom = 1, ni_in_cell
|
||||
i_spec = species_glob(i_atom)
|
||||
do i_atom = 1, n_atoms_pDOS
|
||||
i_spec = species_glob(pDOS_atom_index(i_atom))
|
||||
if(flag_l_resolved .and. flag_lm_resolved) then
|
||||
write(filename,'("Atom",I0.7,"DOS_lm.dat")') i_atom
|
||||
write(filename,'("Atom",I0.7,"DOS_lm.dat")') pDOS_atom_index(i_atom)
|
||||
else if(flag_l_resolved) then
|
||||
write(filename,'("Atom",I0.7,"DOS_l.dat")') i_atom
|
||||
write(filename,'("Atom",I0.7,"DOS_l.dat")') pDOS_atom_index(i_atom)
|
||||
else
|
||||
write(filename,'("Atom",I0.7,"DOS.dat")') i_atom
|
||||
write(filename,'("Atom",I0.7,"DOS.dat")') pDOS_atom_index(i_atom)
|
||||
end if
|
||||
open(unit=17, file=filename)
|
||||
do i_spin = 1, nspin
|
||||
|
|
@ -685,6 +755,109 @@ contains
|
|||
return
|
||||
end subroutine process_pdos
|
||||
|
||||
subroutine process_band_structure
|
||||
|
||||
use datatypes
|
||||
use numbers, ONLY: zero, RD_ERR, twopi, half, one, two, four, six
|
||||
use local, ONLY: eigenvalues, n_bands_total, nkp, wtk, efermi, flag_total_iDOS, &
|
||||
flag_procwf_range_Ef, kx, ky, kz, flag_proc_band_str
|
||||
use read, ONLY: read_eigenvalues, read_psi_coeffs
|
||||
use global_module, ONLY: nspin, n_DOS, E_DOS_min, E_DOS_max, sigma_DOS
|
||||
use units, ONLY: HaToeV
|
||||
|
||||
implicit none
|
||||
|
||||
! Local variables
|
||||
integer :: i_band, i_kp, i_spin, n_DOS_wid, n_band, n_min, n_max, i
|
||||
real(double) :: Ebin, dE_DOS, a, pf_DOS, spin_fac, dE
|
||||
real(double), dimension(nspin) :: total_electrons, iDOS_low
|
||||
real(double), dimension(:,:), allocatable :: total_DOS, iDOS
|
||||
real(double), dimension(:,:), allocatable :: occ
|
||||
|
||||
write(*,fmt='(/2x,"Processing eigenvalues to write band structure")')
|
||||
if(nspin==1) then
|
||||
spin_fac = two
|
||||
else if(nspin==2) then
|
||||
spin_fac = one
|
||||
end if
|
||||
! Read eigenvalues
|
||||
call read_eigenvalues
|
||||
if(abs(E_DOS_min)<RD_ERR) then
|
||||
E_DOS_min = minval(eigenvalues(1,:,:))
|
||||
write(*,fmt='(2x,"Band structure lower limit set automatically: ",f12.5," Ha")') E_DOS_min
|
||||
else
|
||||
write(*,fmt='(2x,"Band structure lower limit set by user: ",f12.5," Ha")') E_DOS_min
|
||||
end if
|
||||
if(abs(E_DOS_max)<RD_ERR) then
|
||||
E_DOS_max = maxval(eigenvalues(n_bands_total,:,:))
|
||||
write(*,fmt='(2x,"Band structure upper limit set automatically: ",f12.5," Ha")') E_DOS_max
|
||||
else
|
||||
write(*,fmt='(2x,"Band structure upper limit set by user: ",f12.5," Ha")') E_DOS_max
|
||||
end if
|
||||
write(*,fmt='(2x,"Writing band structure files")')
|
||||
if(flag_proc_band_str==4) then
|
||||
write(*,fmt='(4x,"X-axis will be k-point index")')
|
||||
else if(flag_proc_band_str==0) then
|
||||
write(*,fmt='(4x,"All k-point coordinates provided")')
|
||||
else if(flag_proc_band_str==1) then
|
||||
write(*,fmt='(4x,"X-axis will be kx")')
|
||||
else if(flag_proc_band_str==2) then
|
||||
write(*,fmt='(4x,"X-axis will be ky")')
|
||||
else if(flag_proc_band_str==3) then
|
||||
write(*,fmt='(4x,"X-axis will be kz")')
|
||||
end if
|
||||
open(unit=17, file="BandStructure.dat")
|
||||
do i_spin = 1, nspin
|
||||
dE = zero
|
||||
if(flag_procwf_range_Ef) dE = -efermi(i_spin)
|
||||
write(17,fmt='("# Spin ",I1)') i_spin
|
||||
write(17,fmt='("# Original Fermi-level: ",f12.5," eV")') HaToeV*efermi(i_spin)
|
||||
write(17,fmt='("# Bands shifted relative to Fermi-level")')
|
||||
do i_band=1,n_bands_total ! All bands
|
||||
if(minval(eigenvalues(i_band, :, i_spin))+dE>=E_DOS_min .and. &
|
||||
maxval(eigenvalues(i_band, :, i_spin))+dE<=E_DOS_max) then
|
||||
write(17,fmt='("# Band ",i6)') i_band
|
||||
do i_kp = 1, nkp
|
||||
if(flag_procwf_range_Ef) then
|
||||
if(flag_proc_band_str==4) then
|
||||
write(17,fmt='(i4,f20.12)') i_kp, &
|
||||
HaToeV*(eigenvalues(i_band, i_kp, i_spin) - efermi(i_spin))
|
||||
else if(flag_proc_band_str==0) then
|
||||
write(17,fmt='(4f20.12)') kx(i_kp), ky(i_kp), kz(i_kp), &
|
||||
HaToeV*(eigenvalues(i_band, i_kp, i_spin) - efermi(i_spin))
|
||||
else if(flag_proc_band_str==1) then
|
||||
write(17,fmt='(2f20.12)') kx(i_kp), &
|
||||
HaToeV*(eigenvalues(i_band, i_kp, i_spin) - efermi(i_spin))
|
||||
else if(flag_proc_band_str==2) then
|
||||
write(17,fmt='(2f20.12)') ky(i_kp), &
|
||||
HaToeV*(eigenvalues(i_band, i_kp, i_spin) - efermi(i_spin))
|
||||
else if(flag_proc_band_str==3) then
|
||||
write(17,fmt='(2f20.12)') kz(i_kp), &
|
||||
HaToeV*(eigenvalues(i_band, i_kp, i_spin) - efermi(i_spin))
|
||||
end if
|
||||
else
|
||||
if(flag_proc_band_str==4) then
|
||||
write(17,fmt='(i4,f20.12)') i_kp, HaToeV*eigenvalues(i_band, i_kp, i_spin)
|
||||
else if(flag_proc_band_str==0) then
|
||||
write(17,fmt='(4f20.12)') kx(i_kp), ky(i_kp), kz(i_kp), &
|
||||
HaToeV*eigenvalues(i_band, i_kp, i_spin)
|
||||
else if(flag_proc_band_str==1) then
|
||||
write(17,fmt='(2f20.12)') kx(i_kp), HaToeV*eigenvalues(i_band, i_kp, i_spin)
|
||||
else if(flag_proc_band_str==2) then
|
||||
write(17,fmt='(2f20.12)') ky(i_kp), HaToeV*eigenvalues(i_band, i_kp, i_spin)
|
||||
else if(flag_proc_band_str==3) then
|
||||
write(17,fmt='(2f20.12)') kz(i_kp), HaToeV*eigenvalues(i_band, i_kp, i_spin)
|
||||
end if
|
||||
end if
|
||||
end do ! i_kp = nkp
|
||||
write(17,fmt='("&")')
|
||||
end if
|
||||
end do
|
||||
end do
|
||||
close(unit=17)
|
||||
return
|
||||
end subroutine process_band_structure
|
||||
|
||||
! Important note
|
||||
!
|
||||
! Formally we have: PAO(\mathbf{r}) = f(r) r^l Y_{lm}(\hat{\mathbf{r}})
|
||||
|
|
|
|||
|
|
@ -11,7 +11,7 @@ contains
|
|||
subroutine read_input
|
||||
|
||||
use global_module, ONLY: flag_assign_blocks, flag_fractional_atomic_coords, nspin, &
|
||||
flag_wf_range_Ef, E_DOS_min, E_DOS_max, sigma_DOS, n_DOS
|
||||
flag_wf_range_Ef, E_DOS_min, E_DOS_max, sigma_DOS, n_DOS, ni_in_cell
|
||||
use local
|
||||
use input_module
|
||||
use numbers
|
||||
|
|
@ -122,6 +122,23 @@ contains
|
|||
i_job = 6
|
||||
else if(leqi(job,'pdo').or.leqi(job,'pro')) then
|
||||
i_job = 7
|
||||
else if(leqi(job,'str').or.leqi(job,'bst')) then
|
||||
i_job = 8
|
||||
! x-axis for band structure plots
|
||||
input_string = fdf_string(1,'Process.BandStrucAxis','n')
|
||||
if(leqi(input_string,'n')) then ! K-point number
|
||||
flag_proc_band_str = 4
|
||||
else if(leqi(input_string,'a')) then ! All k-points
|
||||
flag_proc_band_str = 0
|
||||
else if(leqi(input_string,'x')) then ! kx
|
||||
flag_proc_band_str = 1
|
||||
else if(leqi(input_string,'y')) then ! ky
|
||||
flag_proc_band_str = 2
|
||||
else if(leqi(input_string,'z')) then ! kz
|
||||
flag_proc_band_str = 3
|
||||
else ! All k-points
|
||||
flag_proc_band_str = 4
|
||||
end if
|
||||
end if
|
||||
!
|
||||
charge_stub = fdf_string(80,'Process.ChargeStub','chden')
|
||||
|
|
@ -171,7 +188,7 @@ contains
|
|||
! Energy limits (relative to Ef or absolute)
|
||||
E_wf_min = fdf_double('IO.min_wf_E',zero)
|
||||
E_wf_max = fdf_double('IO.max_wf_E',zero)
|
||||
if(i_job==3 .or. i_job==4 .or. i_job==5) then ! Band-resolved charge or STM
|
||||
if(i_job==3 .or. i_job==4 .or. i_job==5 .or. i_job==8) then ! Band-resolved charge or STM
|
||||
! Read in details of bands output from Conquest
|
||||
n_bands_active=fdf_integer('IO.maxnoWF',0)
|
||||
if(n_bands_active>0) then
|
||||
|
|
@ -194,7 +211,11 @@ contains
|
|||
flag_wf_range = .true.
|
||||
flag_wf_range_Ef = fdf_boolean('IO.WFRangeRelative',.true.)
|
||||
else
|
||||
call cq_abort("No bands specified!")
|
||||
E_wf_min = -BIG
|
||||
E_wf_max = BIG
|
||||
flag_wf_range = .true.
|
||||
flag_wf_range_Ef = fdf_boolean('IO.WFRangeRelative',.true.)
|
||||
write(*,fmt='(2x,"No range specified for bands output; assuming all bands")')
|
||||
end if
|
||||
end if
|
||||
! Now read details of bands to output from processing
|
||||
|
|
@ -234,6 +255,9 @@ contains
|
|||
stop
|
||||
end if
|
||||
flag_by_kpoint = fdf_boolean('Process.outputWF_by_kpoint',.false.)
|
||||
! if output only the real part of WFs (for Gamma-point only)
|
||||
flag_outputWF_real = .false.
|
||||
if (leqi(job,'ban')) flag_outputWF_real = fdf_boolean('Process.outputWF_real',.false.)
|
||||
! DOS
|
||||
! Add flag for window relative to Fermi level
|
||||
E_DOS_min = fdf_double('Process.min_DOS_E',E_wf_min)
|
||||
|
|
@ -252,6 +276,28 @@ contains
|
|||
flag_l_resolved = fdf_boolean('Process.pDOS_l_resolved',.false.)
|
||||
flag_lm_resolved = fdf_boolean('Process.pDOS_lm_resolved',.false.)
|
||||
if(flag_lm_resolved .and. (.not.flag_l_resolved)) flag_l_resolved = .true.
|
||||
! How many atoms?
|
||||
n_atoms_pDOS = fdf_integer('Process.n_atoms_pDOS',0)
|
||||
if(n_atoms_pDOS==0) then ! All atoms
|
||||
n_atoms_pDOS = ni_in_cell
|
||||
allocate(pDOS_atom_index(n_atoms_pDOS))
|
||||
do i=1,ni_in_cell
|
||||
pDOS_atom_index(i)=i
|
||||
end do
|
||||
else
|
||||
allocate(pDOS_atom_index(n_atoms_pDOS))
|
||||
if(fdf_block('pDOS_atoms')) then
|
||||
if(1+block_end-block_start<n_atoms_pDOS) &
|
||||
call cq_abort("Too few atoms in pDOS_atoms: ",&
|
||||
1+block_end-block_start,n_atoms_pDOS)
|
||||
do i=1,n_atoms_pDOS
|
||||
read(unit=input_array(block_start+i-1),fmt=*) pDOS_atom_index(i)
|
||||
end do
|
||||
call fdf_endblock
|
||||
else
|
||||
call cq_abort("Specified n_atoms_pDOS but no pDOS_atoms block")
|
||||
end if
|
||||
end if
|
||||
end if
|
||||
! Now read PS files for atomic information
|
||||
call allocate_species_vars
|
||||
|
|
@ -417,10 +463,10 @@ contains
|
|||
efermi = zero
|
||||
if(nspin==1) then
|
||||
read(17,fmt='(a6,f18.10)') str,efermi(1)
|
||||
write(*,fmt='(4x,"Fermi level: ",f12.5," Ha")') efermi(1)
|
||||
write(*,fmt='(4x,"Fermi level: ",f12.5," Ha (=",f10.3," eV)")') efermi(1), efermi(1)*HaToeV
|
||||
else
|
||||
read(17,fmt='(a6,2f18.10)') str,efermi(1), efermi(2)
|
||||
write(*,fmt='(4x,"Fermi levels: ",2f12.5," Ha")') efermi
|
||||
write(*,fmt='(4x,"Fermi levels: ",2f12.5," Ha (=",2f10.3" eV)")') efermi, efermi*HaToeV
|
||||
end if
|
||||
read(17,*) str
|
||||
! Allocate memory
|
||||
|
|
|
|||
|
|
@ -321,7 +321,7 @@ contains
|
|||
|
||||
use datatypes
|
||||
use numbers
|
||||
use local, ONLY: nkp, wtk, efermi, current, nptsx, nptsy, nptsz, eigenvalues, stm_bias, &
|
||||
use local, ONLY: nkp, wtk, efermi, current, nptsx, nptsy, nptsz, eigenvalues, stm_bias, fermi_offset, &
|
||||
n_bands_total, root_file, grid_z, band_full_to_active
|
||||
use output, ONLY: write_dx_density, write_cube, write_dx_coords
|
||||
use global_module, only : nspin
|
||||
|
|
@ -348,16 +348,20 @@ contains
|
|||
call set_prefac_real(9)
|
||||
allocate(current(nptsx,nptsy,nptsz))
|
||||
allocate(psi(nptsx,nptsy,nptsz))
|
||||
|
||||
if (fermi_offset.ne.zero) write(*,fmt='(2X,"Fermi-offset is ",f10.3," eV")') fermi_offset
|
||||
fermi_offset = fermi_offset/HaToeV ! default of fermi_offset is zero
|
||||
if(stm_bias<0) then
|
||||
Emin = stm_bias/HaToeV + Efermi
|
||||
Emax = Efermi
|
||||
Emin = stm_bias/HaToeV + Efermi + fermi_offset
|
||||
Emax = Efermi + fermi_offset
|
||||
else
|
||||
Emin = Efermi
|
||||
Emax = stm_bias/HaToeV + Efermi
|
||||
Emin = Efermi + fermi_offset
|
||||
Emax = stm_bias/HaToeV + Efermi + fermi_offset
|
||||
end if
|
||||
if(nspin==1) then
|
||||
write(*,fmt='(2x,"Including bands between ",f7.3," and ",f7.3," eV")') Emin(1)*HaToeV, Emax(1)*HaToeV
|
||||
write(*,fmt='(2x,"Including bands between ",f10.3," and ",f10.3," eV")') Emin(1)*HaToeV, Emax(1)*HaToeV
|
||||
end if
|
||||
|
||||
current = zero
|
||||
do ispin=1,nspin
|
||||
do band=1,n_bands_total
|
||||
|
|
|
|||
Loading…
Reference in New Issue