scnc
/
quantum-espresso
mirror of https://gitlab.com/QEF/q-e.git
1
0
Fork 0
Code Issues Packages Projects Releases Wiki Activity

Merge branch 'develop' into 'us_acc' 

# Conflicts:
#   GWW/gww/Makefile
#   GWW/simple_bse/Makefile
#   GWW/simple_ip/Makefile
#   HP/src/hp_solve_linear_system.f90
#   PHonon/PH/solve_e.f90
#   PHonon/PH/solve_linter.f90
#   QEHeat/src/make.depend
This commit is contained in:
giannozz 2021-09-03 19:31:23 +00:00
parent ad8a331aea d9b91704fa
commit 8495cdf4aa
140 changed files with 1832 additions and 1742 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

8
CMakeLists.txt
View File

@ -41,7 +41,6 @@ set(CMAKE_RUNTIME_OUTPUT_DIRECTORY ${qe_BINARY_DIR}/bin
###########################################################
# Build helpers
###########################################################
set(PROJECT_CMAKE ${CMAKE_CURRENT_SOURCE_DIR}/cmake)
set(CMAKE_MODULE_PATH "${CMAKE_CURRENT_SOURCE_DIR}/cmake" ${CMAKE_MODULE_PATH})
include_directories("${CMAKE_CURRENT_SOURCE_DIR}/include")
include(cmake/qeHelpers.cmake)
@ -239,11 +238,11 @@ endif()
# Compiler vendor specific options
############################################################
if(CMAKE_Fortran_COMPILER_ID MATCHES "GNU")
include(${PROJECT_CMAKE}/GNUFortranCompiler.cmake)
include(GNUFortranCompiler)
elseif(CMAKE_Fortran_COMPILER_ID MATCHES "PGI" OR CMAKE_Fortran_COMPILER_ID MATCHES "NVHPC")
include(${PROJECT_CMAKE}/NVFortranCompiler.cmake)
include(NVFortranCompiler)
elseif(CMAKE_Fortran_COMPILER_ID MATCHES "XL")
include(${PROJECT_CMAKE}/IBMFortranCompiler.cmake)
include(IBMFortranCompiler)
endif()
if(QE_ENABLE_STATIC_BUILD)
@ -577,7 +576,6 @@ endif(QE_ENABLE_PROFILE_NVTX)
# Components
###########################################################
add_subdirectory(external)
add_subdirectory(clib)
add_subdirectory(FFTXlib)
add_subdirectory(UtilXlib)
add_subdirectory(Modules)

14
COUPLE/examples/Makefile.gfortran
View File

@ -9,11 +9,7 @@ MPILIBS=-lgfortran -lmpi_mpifh -lmpi
# location of required libraries
PWOBJS = ../src/libqecouple.a ../../PW/src/libpw.a
CPOBJS = ../src/libqecouple.a ../../CPV/src/libcp.a
QEMODS = ../../Modules/libqemod.a ../../KS_Solvers/Davidson/libdavid.a \
../../KS_Solvers/CG/libcg.a ../../FFTXlib/libqefft.a \
../../LAXlib/libqela.a ../../UtilXlib/libutil.a \
../../dft-d3/libdftd3qe.a
LIBOBJS = ../../clib/clib.a ../../iotk/src/libiotk.a
QEMODS=../../KS_Solvers/libks_solvers.a ../../dft-d3/libdftd3qe.a $(BASEMODS)
TLDEPS=couple
@ -21,16 +17,16 @@ TLDEPS=couple
all : tldeps c2pw.x f2pw.x c2cp.x f2cp.x
f2pw.x : f2pw.o $(PWOBJS)
$(LD) $(LDFLAGS) -o $@ $^ $(QEMODS) $(LIBOBJS) $(LIBS)
$(LD) $(LDFLAGS) -o $@ $^ $(QEMODS) $(QELIBS)
c2pw.x : c2pw.o $(PWOBJS)
$(MPICXX) $(LDFLAGS) -o $@ $^ $(MPILIBS) $(QEMODS) $(LIBOBJS) $(LIBS)
$(MPICXX) $(LDFLAGS) -o $@ $^ $(MPILIBS) $(QEMODS) $(QELIBS)
f2cp.x : f2cp.o $(CPOBJS)
$(LD) $(LDFLAGS) -o $@ $^ $(QEMODS) $(LIBOBJS) $(LIBS)
$(LD) $(LDFLAGS) -o $@ $^ $(QEMODS) $(QELIBS)
c2cp.x : c2cp.o $(CPOBJS)
$(MPICXX) $(LDFLAGS) -o $@ $^ $(MPILIBS) $(QEMODS) $(LIBOBJS) $(LIBS)
$(MPICXX) $(LDFLAGS) -o $@ $^ $(MPILIBS) $(QEMODS) $(QELIBS)
%.o: %.cpp
$(MPICXX) -I../include -c $(MPICXXFLAGS) $< -o $@

14
COUPLE/examples/Makefile.ifort
View File

@ -9,11 +9,7 @@ MPILIBS=-Wl,-Bstatic,-lifport,-lifcore,-limf,-Bdynamic -lmpi_mpifh -lmpi
# location of required libraries
PWOBJS = ../src/libqecouple.a ../../PW/src/libpw.a
CPOBJS = ../src/libqecouple.a ../../CPV/src/libcp.a
QEMODS = ../../Modules/libqemod.a ../../KS_Solvers/Davidson/libdavid.a \
../../KS_Solvers/CG/libcg.a ../../FFTXlib/libqefft.a \
../../LAXlib/libqela.a ../../UtilXlib/libutil.a \
../../dft-d3/libdftd3qe.a
LIBOBJS = ../../clib/clib.a ../../iotk/src/libiotk.a
QEMODS=../../KS_Solvers/libks_solvers.a ../../dft-d3/libdftd3qe.a $(BASEMODS)
TLDEPS=couple
@ -21,16 +17,16 @@ TLDEPS=couple
all : tldeps c2pw.x f2pw.x c2cp.x f2cp.x
f2pw.x : f2pw.o $(PWOBJS)
$(LD) $(LDFLAGS) -o $@ $^ $(QEMODS) $(LIBOBJS) $(LIBS)
$(LD) $(LDFLAGS) -o $@ $^ $(QEMODS) $(QELIBS)
c2pw.x : c2pw.o $(PWOBJS)
$(MPICXX) $(LDFLAGS) -o $@ $^ $(MPILIBS) $(QEMODS) $(LIBOBJS) $(LIBS)
$(MPICXX) $(LDFLAGS) -o $@ $^ $(MPILIBS) $(QEMODS) $(QELIBS)
f2cp.x : f2cp.o $(CPOBJS)
$(LD) $(LDFLAGS) -o $@ $^ $(QEMODS) $(LIBOBJS) $(LIBS)
$(LD) $(LDFLAGS) -o $@ $^ $(QEMODS) $(QELIBS)
c2cp.x : c2cp.o $(CPOBJS)
$(MPICXX) $(LDFLAGS) -o $@ $^ $(MPILIBS) $(QEMODS) $(LIBOBJS) $(LIBS)
$(MPICXX) $(LDFLAGS) -o $@ $^ $(MPILIBS) $(QEMODS) $(QELIBS)
%.o: %.cpp
$(MPICXX) -I../include -c $(MPICXXFLAGS) $< -o $@

383
CPV/Doc/user_guide.md
View File

@ -3,42 +3,36 @@
Introduction
============
This guide covers the usage of the `CP` package, version 6.6, a core
This guide covers the usage of the `CP` package, version 6.8, a core
component of the Quantum ESPRESSO distribution. Further documentation,
beyond what is provided in this guide, can be found in the directory
`CPV/Doc/`, containing a copy of this guide.
*Important notice: due to the lack of time and of manpower, this manual
is only partially updated and may contain outdated information.*
This guide assumes that you know the physics that `CP` describes and the
This guide assumes that you know the physics that `CP` describes and the
methods it implements. It also assumes that you have already installed,
or know how to install, Quantum ESPRESSO. If not, please read the
general User's Guide for Quantum ESPRESSO, found in directory `Doc/` two
levels above the one containing this guide; or consult the web site:\
levels above the one containing this guide; or consult the web site:
`http://www.quantum-espresso.org`.
People who want to modify or contribute to `CP` should read the
Developer Manual:\
`Doc/developer_man.pdf`.
People who want to modify or contribute to `CP` should read the
Developer Manual: `https://gitlab.com/QEF/q-e/-/wikis/home`.
`CP` can perform Car-Parrinello molecular dynamics, including
variable-cell dynamics, and free-energy surface calculation at fixed
cell through meta-dynamics, if patched with PLUMED.
The `CP` package is based on the original code written by Roberto Car
and Michele Parrinello. `CP` was developed by Alfredo Pasquarello (EPF
`CP` can perform Car-Parrinello molecular dynamics, including
variable-cell dynamics. The `CP` package is based on the original code
written by Roberto Car
and Michele Parrinello. `CP` was developed by Alfredo Pasquarello (EPF
Lausanne), Kari Laasonen (Oulu), Andrea Trave, Roberto Car (Princeton),
Nicola Marzari (EPF Lausanne), Paolo Giannozzi, and others. FPMD, later
merged with `CP`, was developed by Carlo Cavazzoni, Gerardo Ballabio
(CINECA), Sandro Scandolo (ICTP), Guido Chiarotti, Paolo Focher, and
others. We quote in particular:
merged with `CP`, was developed by Carlo Cavazzoni (Leonardo), Gerardo
Ballabio (CINECA), Sandro Scandolo (ICTP), Guido Chiarotti, Paolo Focher,
and others. We quote in particular:
- Federico Grasselli and Riccardo Bertossa (SISSA) for bug fixes,
extensions to Autopilot;
- Biswajit Santra, Hsin-Yu Ko, Marcus Calegari Andrade (Princeton) for
SCAN functional;
various contribution, notably the SCAN functional;
- Robert DiStasio (Cornell)), Biswajit Santra, and Hsin-Yu Ko for
hybrid functionals with MLWF; (maximally localized Wannier
@ -50,21 +44,21 @@ others. We quote in particular:
- Paolo Umari (Univ. Padua) for finite electric fields and conjugate
gradients;
- Paolo Umari and Ismaila Dabo for ensemble-DFT;
- Paolo Umari and Ismaila Dabo (Penn State) for ensemble-DFT;
- Xiaofei Wang (Princeton) for META-GGA;
- The Autopilot feature was implemented by Targacept, Inc.
This guide has been mostly writen by Gerardo Ballabio and Carlo
Cavazzoni.
The original version of this guide was mostly written by Gerardo Ballabio
and Carlo Cavazzoni.
`CP` is free software, released under the GNU General Public License.\
`CP` is free software, released under the GNU General Public License.\
See `http://www.gnu.org/licenses/old-licenses/gpl-2.0.txt`, or the file
License in the distribution).
`License` in the distribution.
We shall greatly appreciate if scientific work done using the Quantum
ESPRESSO distribution will contain an acknowledgment to the following
ESPRESSO distribution will contain an acknowledgment to the following
references:
> P. Giannozzi, S. Baroni, N. Bonini, M. Calandra, R. Car, C. Cavazzoni,
@ -95,22 +89,22 @@ Users of the GPU-enabled version should also cite the following paper:
> Ferretti, N. Marzari, I. Timrov, A. Urru, S. Baroni, J. Chem. Phys.
> 152, 154105 (2020)
Note the form Quantum ESPRESSO for textual citations of the code. Please
also see package-specific documentation for further recommended
citations. Pseudopotentials should be cited as (for instance)
Note the form `Quantum ESPRESSO` (in small caps) for textual citations
of the code. Please also see other package-specific documentation for
further recommended citations. Pseudopotentials should be cited as
(for instance)
> \[ \] We used the pseudopotentials C.pbe-rrjkus.UPF and O.pbe-vbc.UPF
> from\
> `http://www.quantum-espresso.org`.
> from `http://www.quantum-espresso.org`.
Compilation
===========
`CP` is included in the core Quantum ESPRESSO distribution. Instruction
`CP` is included in the core Quantum ESPRESSO distribution. Instruction
on how to install it can be found in the general documentation (User's
Guide) for Quantum ESPRESSO.
Typing `make cp` from the main Quantum ESPRESSO directory or `make` from
Typing `make cp` from the main Quantum ESPRESSO directory or `make` from
the `CPV/` subdirectory produces the following codes in `CPV/src`:
- `cp.x`: Car-Parrinello Molecular Dynamics code
@ -125,52 +119,44 @@ Symlinks to executable programs will be placed in the `bin/`
subdirectory.
As a final check that compilation was successful, you may want to run
some or all of the tests and examples. Automated tests for `cp.x` are in
some or all of the tests and examples. Automated tests for `cp.x` are in
directory `test-suite/` and can be run via the `Makefile` found there.
Please see the general User's Guide for their setup.
You may take the tests and examples distributed with `CP` as templates
You may take the tests and examples distributed with `CP` as templates
for writing your own input files. Input files for tests are contained in
subdirectories `test-suite/cp_` with file type `*.in1`, `*.in2`, \... .
subdirectories `test-suite/cp_*` with file type `*.in1`, `*.in2`, \... .
Input files for examples are produced, if you run the examples, in the
`results/` subdirectories, with names ending with `.in`.
For general information on parallelism and how to run in parallel
execution, please see the general User's Guide. `CP` currently can take
advantage of both MPI and OpenMP parallelization. The "plane-wave",
"linear-algebra" and "task-group" parallelization levels are
implemented.
execution, please see the general User's Guide. `CP` currently can take
advantage of both MPI and OpenMP parallelization and on GPU acceleration.
The "plane-wave", "linear-algebra" and "task-group" parallelization levels
are implemented.
Input data
==========
Input data for `cp.x` is organized into several namelists, followed by
other fields ("cards") introduced by keywords. The namelists are
other fields ("cards") introduced by keywords. The namelists are:
------------------- ----------------------------------------------------------
&CONTROL: general variables controlling the run
&SYSTEM: structural information on the system under investigation
&ELECTRONS: electronic variables, electron dynamics
&IONS : ionic variables, ionic dynamics
&CELL (optional): variable-cell dynamics
------------------- ----------------------------------------------------------
> &CONTROL: general variables controlling the run\
> &SYSTEM: structural information on the system under investigation\
> &ELECTRONS: electronic variables, electron dynamics\
> &IONS : ionic variables, ionic dynamics\
> &CELL (optional): variable-cell dynamics\
\
The `&CELL` namelist may be omitted for fixed-cell calculations. This
depends on the value of variable `calculation` in namelist &CONTROL.
Most variables in namelists have default values. Only the following
Most variables in namelists have default values. Only he following
variables in &SYSTEM must always be specified:
----------- --------------------- -----------------------------------------------
`ibrav` (integer) Bravais-lattice index
`celldm` (real, dimension 6) crystallographic constants
`nat` (integer) number of atoms in the unit cell
`ntyp` (integer) number of types of atoms in the unit cell
`ecutwfc` (real) kinetic energy cutoff (Ry) for wavefunctions.
----------- --------------------- -----------------------------------------------
\
).
> `ibrav` (integer) Bravais-lattice index\
> `celldm` (real, dimension 6) crystallographic constants\
> `nat` (integer) number of atoms in the unit cell\
> `ntyp` (integer) number of types of atoms in the unit cell\
> `ecutwfc` (real) kinetic energy cutoff (Ry) for wavefunctions
Explanations for the meaning of variables `ibrav` and `celldm`, as well
as on alternative ways to input structural data, are contained in files
@ -178,34 +164,31 @@ as on alternative ways to input structural data, are contained in files
describe a large number of other variables as well. Almost all variables
have default values, which may or may not fit your needs.
Comment lines in namelists can be introduced by a \"!\", exactly as in
fortran code.
After the namelists, you have several fields ("cards") introduced by
keywords with self-explanatory names:
> ATOMIC\_SPECIES\
> ATOMIC\_POSITIONS\
> CELL\_PARAMETERS (optional)\
> OCCUPATIONS (optional)\
> OCCUPATIONS (optional)
The keywords may be followed on the same line by an option. Unknown
fields are ignored. See the files mentioned above for details on the
available "cards".
Comments lines in "cards" can be introduced by either a "!" or a "\#"
character in the first position of a line.
Comment lines in namelists can be introduced by a \"!\", exactly as in
fortran code. Comments lines in "cards" can be introduced by either a "!"
or a "\#" character in the first position of a line.
Data files
----------
The output data files are written in the directory specified by variable
`outdir`, with names specified by variable `prefix` (a string that is
prepended to all file names, whose default value is: `prefix=’pwscf’`).
The `iotk` toolkit is used to write the file in a XML format, whose
definition can be found in the Developer Manual. In order to use the
data directory on a different machine, you need to convert the binary
files to formatted and back, using the `bin/iotk` script.
prepended to all file names, whose default value is `prefix=’cp_$ndw’`,
where `ndw` is an integer specified in input).
In order to use the data on a different machine, you may need to
compile `CP` with HDF5 enabled.
The execution stops if you create a file `prefix.EXIT` either in the
working directory (i.e. where the program is executed), or in the
@ -215,58 +198,13 @@ this procedure is that all files are properly closed, whereas just
killing the process may leave data and output files in an unusable
state.
Format of arrays containing charge density, potential, etc.
-----------------------------------------------------------
The index of arrays used to store functions defined on 3D meshes is
actually a shorthand for three indices, following the FORTRAN convention
(\"leftmost index runs faster\"). An example will explain this better.
Suppose you have a 3D array `psi(nr1x,nr2x,nr3x)`. FORTRAN compilers
store this array sequentially in the computer RAM in the following way:
psi( 1, 1, 1)
psi( 2, 1, 1)
...
psi(nr1x, 1, 1)
psi( 1, 2, 1)
psi( 2, 2, 1)
...
psi(nr1x, 2, 1)
...
...
psi(nr1x,nr2x, 1)
...
psi(nr1x,nr2x,nr3x)
etc
Let `ind` be the position of the `(i,j,k)` element in the above list:
the following relation
ind = i + (j - 1) * nr1x + (k - 1) * nr2x * nr1x
holds. This should clarify the relation between 1D and 3D indexing. In
real space, the `(i,j,k)` point of the FFT grid with dimensions `nr1`
( $`\le`$ `nr1x`), `nr2` ( $`\le`$ `nr2x`), , `nr3` ( $`\le`$ `nr3x`), is
```math
r_{ijk}=\frac{i-1}{nr1} \tau_1 + \frac{j-1}{nr2} \tau_2 + \frac{k-1}{nr3} \tau_3
```
where the $`\tau_i`$ are the basis vectors of the
Bravais lattice. The latter are stored row-wise in the `at` array:
$`\tau_1 =`$ `at(:, 1)`, $`\tau_2 =`$ `at(:, 2)`, $`\tau_3 =`$ `at(:, 3)`.
The distinction between the dimensions of the FFT grid, `(nr1,nr2,nr3)`
and the physical dimensions of the array, `(nr1x,nr2x,nr3x)` is done
only because it is computationally convenient in some cases that the two
sets are not the same. In particular, it is often convenient to have
`nrx1`=`nr1`+1 to reduce memory conflicts.
The format of arrays containing charge density, potential, etc.
is described in the developer manual.
Output files
==========
The `cp.x` code produces many output file, that together build up the trajectory.
The `cp.x` code produces many output files, that together build up the trajectory.
You have a file for the positions, called `prefix.pos`, where `prefix` is defined in
the input file, that is formatted like:
@ -280,35 +218,40 @@ the input file, that is formatted like:
0.42395189282719E+01 0.55766875434652E+01 0.31291744042209E+01
0.45445534106843E+01 0.36049553522533E+01 0.55864387532281E+01
where in the first line there is an header with, in order, the number of the step and
the time in ps of this step. Later you found the positions of all the atoms, in the
same order of the input file (note that this behaviour emerged in v6.6 -- previously
atoms were sorted by type). In this example we have 3 atoms.
The type must be deduced from the input file. After the first 4 lines
you find the same structure for the second step. The units of the position are Bohr's
radius. Note that the atoms coordinates are unwrapped, so it is possible that they go
outside the simulation cell.
where the first line contains the step number and elapsed time, in ps, at this
step; the following lines contain the positions, in Bohr radii, of all the
atoms (3 in this examples), in the same order as in the input file (since v6.6
-- previously, atoms were sorted by type; the type must be deduced from the
input file). The same structure is repeated for the second step and so on.
The printout is made every `iprint` steps (10 in this case, so at step 10, 20,
etc.). Note that the atomic coordinates are not wrapped into the simulation
cell, so it is possible that they lie outside it.
The velocities are written in a similar file named `prefix.vel`, where `prefix` is defined in
the input file, that is formatted like the `.pos` file. The units are the usual Hartree
atomic units (note again that the velocity in the pw code differs by a factor of 2).
The velocities are written in a similar file named `prefix.vel`, where `prefix`
is defined in the input file, that is formatted like the `.pos` file. The units
are the usual Hartree atomic units (note that the velocities in the `pw.x` code
are in _Rydberg_ a.u. and differ by a factor 2).
The `prefix.for` file is formatted like the previous two. Contains the computed forces
and has Hartree atomic units too.
It is written only if `tprnfor = .true.` is set in the input file.
The `prefix.for` file, formatted like the previous two, contains the computed
forces, in Hartree atomic units as well. It is written only if a molecular
dynamics calculation is performed, or if `tprnfor = .true.` is set in input.
The file `prefix.evp` has one line per printed step and contains some thermodynamic data.
The file `prefix.evp` has one line per printed step and contains some
thermodynamical data.
The first line of the file names the columns:
```
# nfi time(ps) ekinc T\_cell(K) Tion(K) etot enthal econs econt Volume Pressure(GPa
# nfi time(ps) ekinc Tcell(K) Tion(K) etot enthal econs econt Volume Pressure(GPa)
```
where:
- `ekinc` $`K_{ELECTRONS}`$, the electron's fake kinetic energy
- `enthal` $`E_{DFT}+PV`$
- `etot` $`E_{DFT}`$ potential energy of the system, the DFT energy
- `econs` $`E_{DFT} + K_{NUCLEI}`$ this is something that is a constant of motion in the limit where the electronic fictitious mass is zero. It has a physical meaning.
- `econt` $`E_{DFT} + K_{IONS} + K_{ELECTRONS}`$ this is a constant of motion of the lagrangian. If the dt is small enough this will be up to a very good precision a constant. It is not a physical quantity, since $`K_{ELECTRONS}`$ has _nothing_ to do with the quantum kinetic energy of the electrons.
- `ekinc` is the electrons fictitious kinetic energy, $`K_{ELECTRONS}`$
- `enthal` is the enthalpy, $`E_{DFT}+PV`$
- `etot` is the DFT (potential) energy of the system, $`E_{DFT}`$
- `econs` is a physically meaningful constant of motion, $`E_{DFT} + K_{NUCLEI}`$,
in the limit of zero electronic fictitious mass
- `econt` is the constant of motion of the lagrangian$`E_{DFT} + K_{IONS} + K_{ELECTRONS}`$ t.
If the time step `dt` is small enough this will be up to a very good precision a constant.
It is not a physical quantity, since $`K_{ELECTRONS}`$ has _nothing_ to do with the quantum
kinetic energy of the electrons.
Using `CP`
@ -317,7 +260,7 @@ Using `CP`
It is important to understand that a CP simulation is a sequence of
different runs, some of them used to \"prepare\" the initial state of
the system, and other performed to collect statistics, or to modify the
state of the system itself, i.e. modify the temperature or the pressure.
state of the system itself, i.e. to modify the temperature or the pressure.
To prepare and run a CP simulation you should first of all define the
system:
@ -393,8 +336,7 @@ An example of input file (Benzene Molecule):
H -2.2 2.2 0.0
H 2.2 2.2 0.0
You can find the description of the input variables in file
`Doc/INPUT_CP.*`.
You can find the description of the input variables in file `Doc/INPUT_CP.*`.
Reaching the electronic ground state
------------------------------------
@ -403,7 +345,7 @@ The first run, when starting from scratch, is always an electronic
minimization, with fixed ions and cell, to bring the electronic system
on the ground state (GS) relative to the starting atomic configuration.
This step is conceptually very similar to self-consistency in a
`pw.x` run.
`pw.x` run.
Sometimes a single run is not enough to reach the GS. In this case, you
need to re-run the electronic minimization stage. Use the input of the
@ -428,14 +370,12 @@ $`< 10^{-5}`$. You could check the value of the fictitious kinetic energy
on the standard output (column EKINC).
Different strategies are available to minimize electrons, but the most
used ones are:
- steepest descent: `electron_dynamics = ’sd’`
- damped dynamics: `electron_dynamics = ’damp’`, `electron_damping` =
a number typically ranging from 0.1 and 0.5
frequently used is _damped dynamics_: `electron_dynamics = ’damp’` and
`electron_damping` = a number typically ranging from 0.1 and 0.5.
See the input description to compute the optimal damping factor.
Steepest descent: `electron_dynamics = ’sd’`, is also available but it
is typicallyslower than damped dynamics and should be used only to
start the minimization.
Relax the system
----------------
@ -860,14 +800,6 @@ ranges between 4 and 7.
All the other parameters have the same meaning in the usual `CP` input,
and they are discussed above.
### Free-energy surface calculations
Once `CP` is patched with `PLUMED` plug-in, it becomes possible to
turn-on most of the PLUMED functionalities running `CP` as:
`./cp.x -plumed` plus the other usual `CP` arguments. The PLUMED input
file has to be located in the specified `outdir` with the fixed name
`plumed.dat`.
### Treatment of USPPs
The cutoff `ecutrho` defines the resolution on the real space FFT mesh
@ -1030,99 +962,62 @@ An example input is listed as following:
O 16.0D0 O_HSCV_PBE-1.0.UPF
H 2.0D0 H_HSCV_PBE-1.0.UPF
Performances
============
Parallel Performances
=====================
`cp.x` can run in principle on any number of processors. The
effectiveness of parallelization is ultimately judged by the "scaling",
i.e. how the time needed to perform a job scales with the number of
processors, and depends upon:
processors. Ideally one would like to have linear scaling, i.e.
$`T \sim T_0/N_p`$ for $`N_p`$ processors, where $`T_0`$ is the estimated
time for serial execution. In addition, one would like to have linear
scaling of the RAM per processor: $`O_N \sim O_0/N_p`$, so that large-memory
systems fit into the RAM of each processor.
- the size and type of the system under study;
We refer to the "Parallelization" section of the general User's Guide for
a description of MPI and OpenMP parallelization paradigms, of the various
MPI parallelization levels, and on how to activate them.
- the judicious choice of the various levels of parallelization
(detailed in
Sec.[\[SubSec:para\]](#SubSec:para){reference-type="ref"
reference="SubSec:para"});
A judicious choice of the various levels of parallelization, together
with the availability of suitable hardware (e.g. fast communications)
is fundamental to reach good performances._VERY IMPORTANT_: For each
system there is an optimal range of number of processors on which to
run the job. A too large number of processors or a bad parallelization
style will yield performance degradation.
- the availability of fast interprocess communications (or lack of
it).
For `CP` with hybrid functionals, see the related section above this one.
For all other cases, the relevant MPI parallelization levels are:
Ideally one would like to have linear scaling, i.e. $`T \sim T_0/N_p`$ for
$`N_p`$ processors, where $`T_0`$ is the estimated time for serial
execution. In addition, one would like to have linear scaling of the RAM
per processor: $`O_N \sim O_0/N_p`$, so that large-memory systems fit into
the RAM of each processor.
- "plane waves" (PW);
- "tasks" (activated by command-line option `-nt N`);
- "linear algebra" (`-nd N`);
- "bands" parallelization (`-nb N`), to be used only in
special cases;
- "images" parallelization (`-ni N`), used only in code `manycp.x`
(see the header of `CPV/src/manycp.f90` for documentation).
As a general rule, image parallelization:
As a rule of thumb:
- start with PW parallelization only (e.g. `mpirun -np N cp.x ...` with
no other parallelization options); the code will scale well unless `N`
exceeds the third FFT dimensions `nr3` and/or `nr3s`.
- To further increase the number of processors, use "task groups",
typically 4 to 8 (e.g. `mpirun -np N cp.x -nt 8 ...`).
- Alternatively, or in addition, you may compile with OpenMP:
`./configure --enable-openmp ...`, then `export OMP_NUM_THREADS=n`
and run on `n` threads (4 to 8 typically).
_Beware conflicts between MPI and OpenMP threads_!
don't do this unless you know what you are doing.
- Finally, the optimal number of processors for \"linear-algebra\"
parallelization can be found by observing the performances of `ortho`
in the final time report for different numbers of processors in the
linear-algebra group (must be a square integer, not larger than the
number of processoris for plane-wave parallelization). Linear-algebra
parallelization distributes `M\times M`$ matrices, with `M` number of
bands, so it may be useful if memory-constrained.
- may give good scaling, but the slowest image will determine the
overall performances ("load balancing" may be a problem);
Note: optimal serial performances are achieved when the data are as much
as possible kept into the cache. As a side effect, PW parallelization may
yield superlinear (better than linear) scaling, thanks to the increase in
serial speed coming from the reduction of data size (making it easier for
the machine to keep data in the cache).
- requires very little communications (suitable for ethernet
communications);
- does not reduce the required memory per processor (unsuitable for
large-memory jobs).
Parallelization on k-points:
- guarantees (almost) linear scaling if the number of k-points is a
multiple of the number of pools;
- requires little communications (suitable for ethernet
communications);
- does not reduce the required memory per processor (unsuitable for
large-memory jobs).
Parallelization on PWs:
- yields good to very good scaling, especially if the number of
processors in a pool is a divisor of $`N_3`$ and $`N_{r3}`$ (the
dimensions along the z-axis of the FFT grids, `nr3` and `nr3s`,
which coincide for NCPPs);
- requires heavy communications (suitable for Gigabit ethernet up to
4, 8 CPUs at most, specialized communication hardware needed for 8
or more processors );
- yields almost linear reduction of memory per processor with the
number of processors in the pool.
A note on scaling: optimal serial performances are achieved when the
data are as much as possible kept into the cache. As a side effect, PW
parallelization may yield superlinear (better than linear) scaling,
thanks to the increase in serial speed coming from the reduction of data
size (making it easier for the machine to keep data in the cache).
VERY IMPORTANT: For each system there is an optimal range of number of
processors on which to run the job. A too large number of processors
will yield performance degradation. If the size of pools is especially
delicate: $`N_p`$ should not exceed $`N_3`$ and $`N_{r3}`$, and should ideally
be no larger than $`1/2\div1/4 N_3`$ and/or $`N_{r3}`$. In order to increase
scalability, it is often convenient to further subdivide a pool of
processors into "task groups". When the number of processors exceeds the
number of FFT planes, data can be redistributed to \"task groups\" so
that each group can process several wavefunctions at the same time.
The optimal number of processors for \"linear-algebra\" parallelization,
taking care of multiplication and diagonalization of $`M\times M`$
matrices, should be determined by observing the performances of
`cdiagh/rdiagh` (`pw.x`) or `ortho` (`cp.x`) for different numbers of
processors in the linear-algebra group (must be a square integer).
Actual parallel performances will also depend on the available software
(MPI libraries) and on the available communication hardware. For PC
clusters, OpenMPI (`http://www.openmpi.org/`) seems to yield better
performances than other implementations (info by Kostantin Kudin). Note
however that you need a decent communication hardware (at least Gigabit
ethernet) in order to have acceptable performances with PW
parallelization. Do not expect good scaling with cheap hardware: PW
calculations are by no means an \"embarrassing parallel\" problem.
Also note that multiprocessor motherboards for Intel Pentium CPUs
typically have just one memory bus for all processors. This dramatically
slows down any code doing massive access to memory (as most codes in the
Quantum ESPRESSO distribution do) that runs on processors of the same
motherboard.

16
CPV/src/Makefile
View File

@ -114,33 +114,33 @@ makov_payne.o
LOBJS = \
entropy.o
QEMODS=../../Modules/libqemod.a ../../upflib/libupf.a ../../FFTXlib/libqefft.a ../../LAXlib/libqela.a ../../UtilXlib/libutil.a ../../XClib/xc_lib.a
QEMODS=$(BASEMODS)
TLDEPS= bindir libs mods
all : tldeps libcp.a manycp.x cp.x wfdd.x cppp.x
manycp.x : manycp.o libcp.a $(QEMODS) $(LIBOBJS)
manycp.x : manycp.o libcp.a $(QEMODS)
$(LD) $(LDFLAGS) -o manycp.x manycp.o \
libcp.a $(QEMODS) $(LIBOBJS) $(QELIBS)
libcp.a $(QEMODS) $(QELIBS)
- ( cd ../../bin ; ln -fs ../CPV/src/manycp.x . )
cp.x : cprstart.o libcp.a $(QEMODS) $(LIBOBJS)
cp.x : cprstart.o libcp.a $(QEMODS)
$(LD) $(LDFLAGS) -o cp.x cprstart.o \
libcp.a $(QEMODS) $(LIBOBJS) $(QELIBS)
libcp.a $(QEMODS) $(QELIBS)
- ( cd ../../bin ; ln -fs ../CPV/src/cp.x . )
libcp.a : $(FOBJS) $(LOBJS)
$(AR) $(ARFLAGS) $@ $?
$(RANLIB) $@
cppp.x : cppp.o $(QEMODS) $(LIBOBJS)
cppp.x : cppp.o $(QEMODS)
$(LD) $(LDFLAGS) -o cppp.x cppp.o \
$(QEMODS) $(LIBOBJS) $(QELIBS)
$(QEMODS) $(QELIBS)
- (cd ../../bin ; ln -fs ../CPV/src/cppp.x . )
wfdd.x : wfdd.o
$(LD) $(LDFLAGS) -o $@ wfdd.o $(QEMODS) $(LIBOBJS) $(QELIBS)
$(LD) $(LDFLAGS) -o $@ wfdd.o $(QEMODS) $(QELIBS)
- ( cd ../../bin ; ln -fs ../CPV/src/$@ . )
tldeps :

35
CPV/src/cg_sub.f90
View File

@ -12,7 +12,8 @@
rhor, rhog, rhos, rhoc, ei1, ei2, ei3, sfac, fion, ema0bg, becdr, &
lambdap, lambda, nlam, vpot, c0, cm, phi, dbec,l_cprestart )
!! please see https://journals.aps.org/rmp/abstract/10.1103/RevModPhys.64.1045
!! please see https://journals.aps.org/prl/abstract/10.1103/PhysRevLett.79.1337 (ensemble DFT)
!! and https://journals.aps.org/rmp/abstract/10.1103/RevModPhys.64.1045 (conjugate gradient)
use kinds, only: dp
use control_flags, only: tpre, iverbosity, tfor, tprnfor
@ -56,7 +57,7 @@
USE cp_main_variables, ONLY : idesc, drhor, drhog
USE mp_global, ONLY: me_image, my_image_id, nbgrp
USE fft_base, ONLY: dffts, dfftp
use wave_gauge, only: project_parallel_gauge_2
!
implicit none
@ -940,33 +941,11 @@
!if required project c0 on previous manifold of occupied states
!NOT IMPLEMENTED YET FOR ENSEMBLE DFT AND NSPIN==2
!NOT IMPLEMENTED FOR US PSEUDOPOTENTIALS
lambda_repl=0.d0
do i = 1, nss
do j = 1, nss
ii = i + istart - 1
jj = j + istart - 1
do ig = 1, ngw
lambda_repl( i, j ) = lambda_repl( i, j ) + &
2.d0 * DBLE( CONJG( c0old( ig, ii ) ) * c0( ig, jj) )
enddo
if( gstart == 2 ) then
lambda_repl( i, j ) = lambda_repl( i, j ) - &
DBLE( CONJG( c0old( 1, ii ) ) * c0( 1, jj ) )
endif
enddo
enddo
CALL mp_sum( lambda_repl, intra_bgrp_comm )
cm(:,:)=c0(:,:)
c0=(0.d0,0.d0)
do i=1,nss
do j=1,nss
c0(1:ngw,i)=c0(1:ngw,i)+lambda_repl(i,j)*cm(1:ngw,j)
enddo
enddo
call project_parallel_gauge_2(c0old, cm, c0, &
nss, ngw, ngw,gstart)
call calbec (nbsp,betae,c0,bec)
CALL gram_bgrp( betae, bec, nkb, c0, ngw )
call calbec(nbsp, betae,c0,bec)

2
CPV/src/chargedensity.f90
View File

@ -1099,7 +1099,7 @@ SUBROUTINE rhov(rhovan,rhog,rhor)
CALL fftx_add_threed2oned_gamma( dfftp, v, rhog(:,isup), rhog(:,isdw) )
!
IF( iverbosity > 1 ) THEN
WRITE( stdout,'(a,2f12.8,/,a,2f12.8)') &
WRITE( stdout,'(a,f12.8,/,a,f12.8)') &
& ' rhov: n_v(g=0) up = ',omega*DBLE (rhog(1,isup)), &
& ' rhov: n_v(g=0) down = ',omega*DBLE(rhog(1,isdw))
END IF

1
CPV/src/make.depend
View File

@ -29,6 +29,7 @@ cg_sub.o : ../../Modules/ions_base.o
cg_sub.o : ../../Modules/kind.o
cg_sub.o : ../../Modules/mp_global.o
cg_sub.o : ../../Modules/recvec.o
cg_sub.o : ../../Modules/wave_gauge.o
cg_sub.o : ../../UtilXlib/mp.o
cg_sub.o : ../../upflib/uspp.o
cg_sub.o : cg.o

2
CPV/src/newd.f90
View File

@ -37,7 +37,7 @@
IMPLICIT NONE
! input
REAL(DP), INTENT(IN) :: rhovan(nhm*(nhm+1)/2,nat,nspin)
REAL(DP) vr(dfftp%nnr,nspin)
REAL(DP), INTENT(IN) :: vr(dfftp%nnr,nspin)
LOGICAL, INTENT(IN) :: tprint
! output
REAL(DP) fion(3,nat)

5
Doc/release-notes
View File

@ -1,8 +1,13 @@
New in development version:
* RMM-DIIS for CPU (S. Nisihara) and GPU (E. de Paoli, P. Delugas)
* DFT-D3: MPI parallelization and GPU acceleration with OPenACC
Fixed in development version:
* Some build problems occurring under special circumstances
* Some PP files were not correctly read since v.6.7
* DFT-D3 with dftd3_version=4 or 6 could produce NaN's in parallel runs
due to missing zero initialization of some work arrays
* Ensemble-DFT in CP ("cg") wasn't working any longer for norm-conserving PPs
Known problems in 6.8 version:
* electron-phonon calculation in the non-colinear/spinorbit case is broken

30
EPW/ZG/src/Makefile
View File

@ -10,42 +10,40 @@ MODFLAGS= $(BASEMOD_FLAGS) \
PHAUXOBJS = ../../../PHonon/PH/libphaux.a
PHOBJS = ../../../PHonon/PH/libph.a
PWOBJS = ../../../PW/src/libpw.a
PWOBJS = ../../../PW/src/libpw.a ../../../KS_Solvers/libks_solvers.a ../../../dft-d3/libdftd3qe.a
LRMODS = ../../../LR_Modules/liblrmod.a
QEMODS = ../../../Modules/libqemod.a ../../../KS_Solvers/libks_solvers.a ../../../upflib/libupf.a ../../../XClib/xc_lib.a \
../../../FFTXlib/libqefft.a ../../../LAXlib/libqela.a ../../../UtilXlib/libutil.a ../../../dft-d3/libdftd3qe.a
LIBOBJS = ../../../clib/clib.a
QEMODS = $(BASEMODS)
all : ZG.x disca.x pp_disca.x pp_spctrlfn.x bands_unfold.x epsilon_Gaus.x
ZG.x : ZG.o $(PHAUXOBJS) $(PHOBJS) $(PWOBJS) $(LRMODS) $(QEMODS) $(LIBOBJS)
ZG.x : ZG.o $(PHAUXOBJS) $(PHOBJS) $(PWOBJS) $(LRMODS) $(QEMODS)
$(LD) $(LDFLAGS) -o $@ \
ZG.o $(PHAUXOBJS) $(PHOBJS) $(LRMODS) $(PWOBJS) $(QEMODS) $(LIBOBJS) $(QELIBS)
ZG.o $(PHAUXOBJS) $(PHOBJS) $(LRMODS) $(PWOBJS) $(QEMODS) $(QELIBS)
- ( cd ../../../bin ; ln -fs ../EPW/ZG/src/$@ . )
disca.x : disca.o $(PHAUXOBJS) $(PHOBJS) $(PWOBJS) $(LRMODS) $(QEMODS) $(LIBOBJS)
disca.x : disca.o $(PHAUXOBJS) $(PHOBJS) $(PWOBJS) $(LRMODS) $(QEMODS)
$(LD) $(LDFLAGS) -o $@ \
disca.o $(PHAUXOBJS) $(PHOBJS) $(LRMODS) $(PWOBJS) $(QEMODS) $(LIBOBJS) $(QELIBS)
disca.o $(PHAUXOBJS) $(PHOBJS) $(LRMODS) $(PWOBJS) $(QEMODS) $(QELIBS)
- ( cd ../../../bin ; ln -fs ../EPW/ZG/src/$@ . )
pp_disca.x : pp_disca.o $(PWOBJS) $(QEMODS) $(LIBOBJS)
pp_disca.x : pp_disca.o $(PWOBJS) $(QEMODS)
$(LD) $(LDFLAGS) -o $@ \
pp_disca.o $(PHAUXOBJS) $(PHOBJS) $(PWOBJS) $(QEMODS) $(LIBOBJS) $(QELIBS)
pp_disca.o $(PHAUXOBJS) $(PHOBJS) $(PWOBJS) $(QEMODS) $(QELIBS)
- ( cd ../../../bin ; ln -fs ../EPW/ZG/src/$@ . )
bands_unfold.x : bands_unfold.o $(PWOBJS) $(QEMODS) $(LIBOBJS)
bands_unfold.x : bands_unfold.o $(PWOBJS) $(QEMODS)
$(LD) $(LDFLAGS) -o $@ \
bands_unfold.o ../../../PP/src/libpp.a $(PWOBJS) $(QEMODS) $(LIBOBJS) $(QELIBS)
bands_unfold.o ../../../PP/src/libpp.a $(PWOBJS) $(QEMODS) $(QELIBS)
- ( cd ../../../bin ; ln -fs ../EPW/ZG/src/$@ . )
pp_spctrlfn.x : pp_spctrlfn.o $(PWOBJS) $(QEMODS) $(LIBOBJS)
pp_spctrlfn.x : pp_spctrlfn.o $(PWOBJS) $(QEMODS)
$(LD) $(LDFLAGS) -o $@ \
pp_spctrlfn.o $(PHAUXOBJS) $(PHOBJS) $(PWOBJS) $(QEMODS) $(LIBOBJS) $(QELIBS)
pp_spctrlfn.o $(PHAUXOBJS) $(PHOBJS) $(PWOBJS) $(QEMODS) $(QELIBS)
- ( cd ../../../bin ; ln -fs ../EPW/ZG/src/$@ . )
epsilon_Gaus.x : epsilon_Gaus.o $(PWOBJS) $(QEMODS) $(LIBOBJS)
epsilon_Gaus.x : epsilon_Gaus.o $(PWOBJS) $(QEMODS)
$(LD) $(LDFLAGS) -o $@ \
epsilon_Gaus.o ../../../PP/src/libpp.a $(PWOBJS) $(QEMODS) $(LIBOBJS) $(QELIBS)
epsilon_Gaus.o ../../../PP/src/libpp.a $(PWOBJS) $(QEMODS) $(QELIBS)
- ( cd ../../../bin ; ln -fs ../EPW/ZG/src/$@ . )
clean :

11
EPW/src/Makefile
View File

@ -78,13 +78,10 @@ EPWOBJS += ephblochkq.o wfc_elec.o test_tools.o
PHOBJS = ../../PHonon/PH/libph.a
PWOBJS = ../../PW/src/libpw.a
W90LIB = ../../wannier90-3.1.0/libwannier.a
LRMODS = ../../LR_Modules/liblrmod.a
PWOBJS = ../../PW/src/libpw.a
QEMODS = ../../Modules/libqemod.a ../../KS_Solvers/libks_solvers.a ../../XClib/xc_lib.a \
../../upflib/libupf.a ../../FFTXlib/libqefft.a ../../dft-d3/libdftd3qe.a
LIBOBJS =../../LAXlib/libqela.a ../../UtilXlib/libutil.a ../../clib/clib.a
PWOBJS = ../../PW/src/libpw.a ../../KS_Solvers/libks_solvers.a ../../dft-d3/libdftd3qe.a
QEMODS = $(BASEMODS)
TLDEPS= pw ph pp
@ -94,9 +91,9 @@ libepw.a : $(EPWOBJS)
$(AR) $(ARFLAGS) $@ $?
$(RANLIB) $@
epw.x : epw.o libepw.a $(PHOBJS) $(LRMODS) $(PWOBJS) $(W90LIB) $(QEMODS) $(LIBOBJS)
epw.x : epw.o libepw.a $(PHOBJS) $(LRMODS) $(PWOBJS) $(W90LIB) $(QEMODS)
$(LD) $(LDFLAGS) -o $@ \
epw.o libepw.a $(PHOBJS) $(LRMODS) $(W90LIB) $(PWOBJS) $(QEMODS) $(LIBOBJS) $(QELIBS)
epw.o libepw.a $(PHOBJS) $(LRMODS) $(W90LIB) $(PWOBJS) $(QEMODS) $(QELIBS)
- ( cd ../bin ; ln -fs ../src/epw.x . )
pw :

4
EPW/src/close_epw.f90
View File

@ -130,8 +130,8 @@
SUBROUTINE close_final
!------------------------------------------------------------------
!
USE units_lr, ONLY : iuwfc
USE units_ph, ONLY : iudwf, iudrho
USE units_lr, ONLY : iuwfc, iudwf
USE units_ph, ONLY : iudrho
USE phcom, ONLY : fildrho
USE mp_global, ONLY : me_pool,root_pool
USE io_var, ONLY : iunepmatwe

106
FFTXlib/fft_scalar.FFTW3.f90
View File

@ -160,16 +160,22 @@
SUBROUTINE init_plan()
implicit none
!
COMPLEX(DP), ALLOCATABLE :: c_test(:)
!
ALLOCATE(c_test, mold=c)
!
IF( C_ASSOCIATED(fw_planz( icurrent)) ) CALL fftw_destroy_plan( fw_planz( icurrent) )
IF( C_ASSOCIATED(bw_planz( icurrent)) ) CALL fftw_destroy_plan( bw_planz( icurrent) )
idir = -1
fw_planz(icurrent) = fftw_plan_many_dft(1, (/nz/), nsl, c, &
(/SIZE(c)/), 1, ldz, cout, (/SIZE(cout)/), 1, ldz, idir, FFTW_ESTIMATE)
fw_planz(icurrent) = fftw_plan_many_dft(1, (/nz/), nsl, c_test, &
(/SIZE(c)/), 1, ldz, cout, (/SIZE(cout)/), 1, ldz, idir, FFTW_MEASURE)
idir = 1
bw_planz(icurrent) = fftw_plan_many_dft(1, (/nz/), nsl, c, &
(/SIZE(c)/), 1, ldz, cout, (/SIZE(cout)/), 1, ldz, idir, FFTW_ESTIMATE)
bw_planz(icurrent) = fftw_plan_many_dft(1, (/nz/), nsl, c_test, &
(/SIZE(c)/), 1, ldz, cout, (/SIZE(cout)/), 1, ldz, idir, FFTW_MEASURE)
!
DEALLOCATE(c_test)
!
zdims(1,icurrent) = nz; zdims(2,icurrent) = nsl; zdims(3,icurrent) = ldz;
ip = icurrent
icurrent = MOD( icurrent, ndims ) + 1
@ -309,42 +315,47 @@
SUBROUTINE init_plan()
implicit none
COMPLEX(DP), ALLOCATABLE :: f_test(:)
!
ALLOCATE(f_test,mold=r)
!
IF ( ldx /= nx .OR. ldy /= ny ) THEN
IF( C_ASSOCIATED(fw_plan(2,icurrent)) ) CALL fftw_destroy_plan( fw_plan(2,icurrent) )
IF( C_ASSOCIATED(bw_plan(2,icurrent)) ) CALL fftw_destroy_plan( bw_plan(2,icurrent) )
idir = -1
fw_plan(2,icurrent) = fftw_plan_many_dft(1, (/ny/), 1, r(1:), &
(/ldx*ldy/), ldx, 1, r(1:), (/ldx*ldy/), ldx, 1, idir, &
FFTW_ESTIMATE)
fw_plan(2,icurrent) = fftw_plan_many_dft(1, (/ny/), 1, f_test(1:), &
(/ldx*ldy/), ldx, 1, f_test(1:), (/ldx*ldy/), ldx, 1, idir, &
FFTW_MEASURE)
idir = 1
bw_plan(2,icurrent) = fftw_plan_many_dft(1, (/ny/), 1, r(1:), &
(/ldx*ldy/), ldx, 1, r(1:), (/ldx*ldy/), ldx, 1, idir, &
FFTW_ESTIMATE)
bw_plan(2,icurrent) = fftw_plan_many_dft(1, (/ny/), 1, f_test(1:), &
(/ldx*ldy/), ldx, 1, f_test(1:), (/ldx*ldy/), ldx, 1, idir, &
FFTW_MEASURE)
IF( C_ASSOCIATED(fw_plan(1,icurrent)) ) CALL fftw_destroy_plan( fw_plan(1,icurrent) )
IF( C_ASSOCIATED(bw_plan(1,icurrent)) ) CALL fftw_destroy_plan( bw_plan(1,icurrent) )
idir = -1
fw_plan(1,icurrent) = fftw_plan_many_dft(1, (/nx/), ny, r(1:), &
(/ldx*ldy/), 1, ldx, r(1:), (/ldx*ldy/), 1, ldx, idir, &
FFTW_ESTIMATE)
fw_plan(1,icurrent) = fftw_plan_many_dft(1, (/nx/), ny, f_test(1:), &
(/ldx*ldy/), 1, ldx, f_test(1:), (/ldx*ldy/), 1, ldx, idir, &
FFTW_MEASURE)
idir = 1
bw_plan(1,icurrent) = fftw_plan_many_dft(1, (/nx/), ny, r(1:), &
(/ldx*ldy/), 1, ldx, r(1:), (/ldx*ldy/), 1, ldx, idir, &
FFTW_ESTIMATE)
bw_plan(1,icurrent) = fftw_plan_many_dft(1, (/nx/), ny, f_test(1:), &
(/ldx*ldy/), 1, ldx, f_test(1:), (/ldx*ldy/), 1, ldx, idir, &
FFTW_MEASURE)
ELSE
IF( C_ASSOCIATED(fw_plan( 1, icurrent)) ) CALL fftw_destroy_plan( fw_plan( 1, icurrent) )
IF( C_ASSOCIATED(bw_plan( 1, icurrent)) ) CALL fftw_destroy_plan( bw_plan( 1, icurrent) )
idir = -1
fw_plan(1, icurrent) = fftw_plan_many_dft(2, (/ny, nx/), nzl,&
r(1:), (/ny, nx/), 1, nx*ny, r(1:), (/ny, nx/), 1, nx*ny, idir,&
FFTW_ESTIMATE)
f_test(1:), (/ldy, ldx/), 1, ldx*ldy, f_test(1:), (/ldy, ldx/), 1, ldx*ldy, idir,&
FFTW_MEASURE)
idir = 1
bw_plan(1, icurrent) = fftw_plan_many_dft(2, (/ny, nx/), nzl,&
r(1:), (/ny, nx/), 1, nx*ny, r(1:), (/ny, nx/), 1, nx*ny, idir,&
FFTW_ESTIMATE)
f_test(1:), (/ldy, ldx/), 1, ldx*ldy, f_test(1:), (/ldy, ldx/), 1, ldx*ldy, idir,&
FFTW_MEASURE)
END IF
!
DEALLOCATE(f_test)
!
dims(1,icurrent) = ny; dims(2,icurrent) = ldx;
dims(3,icurrent) = nx; dims(4,icurrent) = nzl;
ip = icurrent
@ -450,15 +461,21 @@
SUBROUTINE init_plan()
implicit none
COMPLEX(DP), ALLOCATABLE :: f_test(:)
IF ( nx /= ldx .or. ny /= ldy .or. nz /= ldz ) &
call fftx_error__('cfft3','not implemented',3)
IF( C_ASSOCIATED(fw_plan(icurrent)) ) CALL fftw_destroy_plan( fw_plan(icurrent) )
IF( C_ASSOCIATED(bw_plan(icurrent)) ) CALL fftw_destroy_plan( bw_plan(icurrent) )
!
ALLOCATE(f_test,mold=f)
!
idir = -1
fw_plan(icurrent) = fftw_plan_dft_3d(nz, ny, nx, f(1:), f(1:), idir, FFTW_ESTIMATE)
fw_plan(icurrent) = fftw_plan_dft_3d(nz, ny, nx, f_test(1:), f_test(1:), idir, FFTW_MEASURE)
idir = 1
bw_plan(icurrent) = fftw_plan_dft_3d(nz, ny, nx, f(1:), f(1:), idir, FFTW_ESTIMATE)
bw_plan(icurrent) = fftw_plan_dft_3d(nz, ny, nx, f_test(1:), f_test(1:), idir, FFTW_MEASURE)
!
DEALLOCATE(f_test)
!
dims(1,icurrent) = nx; dims(2,icurrent) = ny; dims(3,icurrent) = nz
ip = icurrent
icurrent = MOD( icurrent, ndims ) + 1
@ -630,7 +647,9 @@ SUBROUTINE cfft3ds (f, nx, ny, nz, ldx, ldy, ldz, howmany, isign, &
SUBROUTINE init_plan()
implicit none
!
COMPLEX(DP), ALLOCATABLE :: f_test(:)
!
IF( C_ASSOCIATED(fw_plan( 1, icurrent)) ) &
CALL fftw_destroy_plan( fw_plan( 1, icurrent) )
IF( C_ASSOCIATED(bw_plan( 1, icurrent)) ) &
@ -643,25 +662,30 @@ SUBROUTINE cfft3ds (f, nx, ny, nz, ldx, ldy, ldz, howmany, isign, &
CALL fftw_destroy_plan( fw_plan( 3, icurrent) )
IF( C_ASSOCIATED(bw_plan( 3, icurrent)) ) &
CALL fftw_destroy_plan( bw_plan( 3, icurrent) )
!
ALLOCATE(f_test, mold=f)
!
idir = -1
fw_plan(1, icurrent) = fftw_plan_many_dft(1, (/nx/), ny*nz, f(1:), (/ldz, ldy, ldx/), 1, ldx, &
f(1:), (/ldz, ldy, ldx/), 1, ldx, idir, FFTW_ESTIMATE)
fw_plan(1, icurrent) = fftw_plan_many_dft(1, (/nx/), ny*nz, f_test(1:), (/ldz, ldy, ldx/), 1, ldx, &
f_test(1:), (/ldz, ldy, ldx/), 1, ldx, idir, FFTW_MEASURE)
idir = 1
bw_plan(1, icurrent) = fftw_plan_many_dft(1, (/nx/), ny*nz, f(1:), (/ldz, ldy, ldx/), 1, ldx, &
f(1:), (/ldz, ldy, ldx/), 1, ldx, idir, FFTW_ESTIMATE)
bw_plan(1, icurrent) = fftw_plan_many_dft(1, (/nx/), ny*nz, f_test(1:), (/ldz, ldy, ldx/), 1, ldx, &
f_test(1:), (/ldz, ldy, ldx/), 1, ldx, idir, FFTW_MEASURE)
idir = -1
fw_plan(2, icurrent) = fftw_plan_many_dft(1, (/ny/), nz, f(1:), (/ldz, ldy, ldx/), ldx, ldx*ldy, &
f(1:), (/ldz, ldy, ldx/), ldx, ldx*ldy, idir, FFTW_ESTIMATE)
fw_plan(2, icurrent) = fftw_plan_many_dft(1, (/ny/), nz, f_test(1:), (/ldz, ldy, ldx/), ldx, ldx*ldy, &
f_test(1:), (/ldz, ldy, ldx/), ldx, ldx*ldy, idir, FFTW_MEASURE)
idir = 1
bw_plan(2, icurrent) = fftw_plan_many_dft(1, (/ny/), nz, f(1:), (/ldz, ldy, ldx/), ldx, ldx*ldy, &
f(1:), (/ldz, ldy, ldx/), ldx, ldx*ldy, idir, FFTW_ESTIMATE)
bw_plan(2, icurrent) = fftw_plan_many_dft(1, (/ny/), nz, f_test(1:), (/ldz, ldy, ldx/), ldx, ldx*ldy, &
f_test(1:), (/ldz, ldy, ldx/), ldx, ldx*ldy, idir, FFTW_MEASURE)
idir = -1
fw_plan(3, icurrent) = fftw_plan_many_dft(1, (/nz/), 1, f(1:), (/ldz, ldy, ldx/), ldx*ldy, 1, &
f(1:), (/ldz, ldy, ldx/), ldx*ldy, 1, idir, FFTW_ESTIMATE)
fw_plan(3, icurrent) = fftw_plan_many_dft(1, (/nz/), 1, f_test(1:), (/ldz, ldy, ldx/), ldx*ldy, 1, &
f_test(1:), (/ldz, ldy, ldx/), ldx*ldy, 1, idir, FFTW_MEASURE)
idir = 1
bw_plan(3, icurrent) = fftw_plan_many_dft(1, (/nz/), 1, f(1:), (/ldz, ldy, ldx/), ldx*ldy, 1, &
f(1:), (/ldz, ldy, ldx/), ldx*ldy, 1, idir, FFTW_ESTIMATE)
bw_plan(3, icurrent) = fftw_plan_many_dft(1, (/nz/), 1, f_test(1:), (/ldz, ldy, ldx/), ldx*ldy, 1, &
f_test(1:), (/ldz, ldy, ldx/), ldx*ldy, 1, idir, FFTW_MEASURE)
!
DEALLOCATE(f_test)
!
dims(1,icurrent) = nx; dims(2,icurrent) = ny; dims(3,icurrent) = nz
ip = icurrent
icurrent = MOD( icurrent, ndims ) + 1

105
FFTXlib/fft_test.f90
View File

@ -67,7 +67,9 @@ program test
!!
!!-ntg Number of task groups
!!
!!-gamma Enables gamma point trick. Should be about 2 times faster.
!!-gamma Enables gamma point trick. Should be about 2 times faster
!!
!!-pd If .true. uses pencil decomposition, otherwise uses slab decomposition
!!
!!-av1 x y z First lattice vector, in atomic units. N.B.: when using -av1, -alat is ignored!
!!
@ -90,8 +92,19 @@ program test
USE fft_helper_subroutines
USE fft_interfaces, ONLY:fwfft, invfft
USE timers
!
IMPLICIT NONE
!
ENUM, BIND(C)
ENUMERATOR :: STAMP_BEGIN = 1
ENUMERATOR :: STAMP_PSI = 2
ENUMERATOR :: STAMP_INVFFT = 3
ENUMERATOR :: STAMP_VLOC = 4
ENUMERATOR :: STAMP_FWFFT = 5
ENUMERATOR :: STAMP_HPSI = 6
END ENUM
INTEGER, PARAMETER :: NUM_STAMPS = 6
!
TYPE(fft_type_descriptor) :: dfftp, dffts, dfft3d
!
TYPE(sticks_map) :: smap
@ -121,17 +134,21 @@ program test
!! cut-off for the wave-function
REAL*8 :: tpiba, alat, alat_in
!! lattice parameters
REAL*8 :: time(100)
REAL*8 :: my_time(100)
REAL*8 :: time_min(100)
REAL*8 :: time_max(100)
REAL*8 :: time_avg(100)
REAL*8 :: time(NUM_STAMPS)
REAL*8 :: my_time(NUM_STAMPS)
REAL*8 :: time_min(NUM_STAMPS)
REAL*8 :: time_max(NUM_STAMPS)
REAL*8 :: time_avg(NUM_STAMPS)
REAL*8 :: wall
REAL*8 :: wall_avg
!
LOGICAL :: gamma_only = .false.
LOGICAL :: use_tg
!! if calculations require only gamma point
LOGICAL :: use_tg
!! if calculations use task group
LOGICAL :: use_pd = .false.
!! if calculations use pencil decomposition
LOGICAL :: lpara
REAL*8 :: at(3, 3), bg(3, 3)
REAL(DP), PARAMETER :: pi = 4.0_DP * atan(1.0_DP)
!
@ -247,6 +264,10 @@ program test
CALL get_command_argument(i + 1, arg)
READ (arg, *) gamma_only
END IF
IF (TRIM(arg) == '-pd') THEN
CALL get_command_argument(i + 1, arg)
READ (arg, *) use_pd
END IF
IF ((TRIM(arg) == '-howmany').or.(TRIM(arg) == '-nh')) THEN
CALL get_command_argument(i + 1, arg)
READ (arg, *) many_fft
@ -357,6 +378,7 @@ program test
write (*, *) 'Num Task Group = ', ntgs
write (*, *) 'Num Many FFTs = ', many_fft
write (*, *) 'Gamma trick = ', gamma_only
write (*, *) 'Pencil decomp = ', use_pd
end if
!
nx = 2*int(sqrt(gcutm)*sqrt(at(1, 1)**2 + at(2, 1)**2 + at(3, 1)**2)) + 1
@ -371,12 +393,16 @@ program test
IF (gamma_only) incr = 2
dffts%has_task_groups = (ntgs > 1)
use_tg = dffts%has_task_groups
lpara = (npes > 1)
!
dffts%rho_clock_label='ffts' ; dffts%wave_clock_label='fftw'
CALL fft_type_init(dffts, smap, "wave", gamma_only, .true., comm, at, bg, gkcut, gcutms/gkcut, nyfft=ntgs, nmany=many_fft)
dffts%rho_clock_label='ffts'
dffts%wave_clock_label='fftw'
!
CALL fft_type_init(dffts, smap, "wave", gamma_only, lpara, comm, at, bg, gkcut, gcutms/gkcut, &
nyfft=ntgs, nmany=many_fft, use_pd=use_pd)
dfftp%rho_clock_label='fft'
CALL fft_type_init(dfftp, smap, "rho", gamma_only, .true., comm, at, bg, gcutm, 4.d0, nyfft=ntgs, nmany=many_fft)
CALL fft_type_init(dfftp, smap, "rho", gamma_only, lpara, comm, at, bg, gcutm, 4.d0, nyfft=ntgs, &
nmany=many_fft, use_pd=use_pd)
!
CALL fft_base_info(mype == 0, dffts, dfftp)
if (mype == 0) then
@ -431,7 +457,7 @@ program test
g, gg, mill, ig_l2g, gstart, .TRUE. )
ELSE
CALL ggen( dfftp, gamma_only, at, bg, gcutm, ngm_g, ngm, &
g, gg, mill, ig_l2g, gstart, .FALSE. )
g, gg, mill, ig_l2g, gstart, .FALSE. )
END IF
CALL ggens( dffts, gamma_only, at, g, gg, mill, gcutms, ngms )
!
@ -500,13 +526,13 @@ program test
IF (use_tg) THEN
DO ib = 1, nbnd, incr
!
time(1) = mpi_wall_time()
time(STAMP_BEGIN) = mpi_wall_time()
!
call prepare_psi_tg(ib, nbnd, ngms, psi, tg_psic, dffts, gamma_only)
time(2) = mpi_wall_time()
time(STAMP_PSI) = mpi_wall_time()
!
CALL invfft('tgWave', tg_psic, dffts);
time(3) = mpi_wall_time()
CALL invfft('tgWave', tg_psic, dffts)
time(STAMP_INVFFT) = mpi_wall_time()
!
CALL tg_get_group_nr3(dffts, right_nr3)
!
@ -514,13 +540,13 @@ program test
tg_psic(j) = tg_psic(j)*tg_v(j)
ENDDO
!
time(4) = mpi_wall_time()
time(STAMP_VLOC) = mpi_wall_time()
!
CALL fwfft('tgWave', tg_psic, dffts);
time(5) = mpi_wall_time()
CALL fwfft('tgWave', tg_psic, dffts)
time(STAMP_FWFFT) = mpi_wall_time()
!
CALL accumulate_hpsi_tg(ib, nbnd, ngms, hpsi, tg_psic, dffts, gamma_only)
time(6) = mpi_wall_time()
time(STAMP_HPSI) = mpi_wall_time()
!
DO i = 2, 6
my_time(i) = my_time(i) + (time(i) - time(i - 1))
@ -531,6 +557,8 @@ program test
ENDDO
ELSEIF (many_fft > 1) THEN
DO ib = 1, nbnd, many_fft
!
time(STAMP_BEGIN) = mpi_wall_time()
!
group_size = MIN(many_fft, nbnd - (ib -1))
!
@ -538,24 +566,24 @@ program test
!call prepare_psi(ib, nbnd, ngms, psi, psic(1+(k-1)*dffts%nnr:), dffts, gamma_only)
call prepare_psi(ib, nbnd, ngms, psi, psic, dffts, gamma_only)
ENDDO
time(2) = mpi_wall_time()
time(STAMP_PSI) = mpi_wall_time()
!
CALL invfft('Wave', psic, dffts, howmany=group_size)
time(3) = mpi_wall_time()
time(STAMP_INVFFT) = mpi_wall_time()
!
DO j = 1, dffts%nnr
psic(j) = psic(j)*v(j)
ENDDO
time(4) = mpi_wall_time()
time(STAMP_VLOC) = mpi_wall_time()
!
CALL fwfft('Wave', psic, dffts, howmany=group_size)
time(5) = mpi_wall_time()
time(STAMP_FWFFT) = mpi_wall_time()
!
DO k=0, group_size - 1
!CALL accumulate_hpsi(ib, nbnd, ngms, hpsi, psic(1+(k-1)*dffts%nnr:), dffts, gamma_only)
CALL accumulate_hpsi(ib, nbnd, ngms, hpsi, psic, dffts, gamma_only)
ENDDO
time(6) = mpi_wall_time()
time(STAMP_HPSI) = mpi_wall_time()
!
DO i = 2, 6
my_time(i) = my_time(i) + (time(i) - time(i - 1))
@ -567,21 +595,24 @@ program test
ELSE
DO ib = 1, nbnd, incr
!
call prepare_psi(ib, nbnd, ngms, psi, psic, dffts, gamma_only)
time(2) = mpi_wall_time()
time(STAMP_BEGIN) = mpi_wall_time()
!
CALL invfft('Wave', psic, dffts); time(3) = mpi_wall_time()
call prepare_psi(ib, nbnd, ngms, psi, psic, dffts, gamma_only)
time(STAMP_PSI) = mpi_wall_time()
!
CALL invfft('Wave', psic, dffts)
time(STAMP_INVFFT) = mpi_wall_time()
!
DO j = 1, dffts%nnr
psic(j) = psic(j)*v(j)
ENDDO
time(4) = mpi_wall_time()
time(STAMP_VLOC) = mpi_wall_time()
!
CALL fwfft('Wave', psic, dffts);
time(5) = mpi_wall_time()
CALL fwfft('Wave', psic, dffts)
time(STAMP_FWFFT) = mpi_wall_time()
!
CALL accumulate_hpsi(ib, nbnd, ngms, hpsi, psic, dffts, gamma_only)
time(6) = mpi_wall_time()
time(STAMP_HPSI) = mpi_wall_time()
!
DO i = 2, 6
my_time(i) = my_time(i) + (time(i) - time(i - 1))
@ -639,11 +670,11 @@ program test
write(*,100)
write(*,1)
write(*,100)
write(*,2) time_min(2), time_max(2), time_avg(2)
write(*,3) time_min(3), time_max(3), time_avg(3)
write(*,4) time_min(4), time_max(4), time_avg(4)
write(*,5) time_min(5), time_max(5), time_avg(5)
write(*,6) time_min(6), time_max(6), time_avg(6)
write(*,2) time_min(STAMP_PSI), time_max(STAMP_PSI), time_avg(STAMP_PSI)
write(*,3) time_min(STAMP_INVFFT), time_max(STAMP_INVFFT), time_avg(STAMP_INVFFT)
write(*,4) time_min(STAMP_VLOC), time_max(STAMP_VLOC), time_avg(STAMP_VLOC)
write(*,5) time_min(STAMP_FWFFT), time_max(STAMP_FWFFT), time_avg(STAMP_FWFFT)
write(*,6) time_min(STAMP_HPSI), time_max(STAMP_HPSI), time_avg(STAMP_HPSI)
write(*,7) wall
write(*,100)

8
GWW/bse/Makefile
View File

@ -45,10 +45,8 @@ qpcorrections.o
#zvscal_test.o\
QEMODS = ../../Modules/libqemod.a ../../FFTXlib/libqefft.a \
../../upflib/libupf.a ../../KS_Solvers/libks_solvers.a \
../../LAXlib/libqela.a ../../UtilXlib/libutil.a ../../XClib/xc_lib.a
PWOBJS = ../../PW/src/libpw.a
QEMODS = $(BASEMODS)
PWOBJS = ../../PW/src/libpw.a ../../KS_Solvers/libks_solvers.a
PW4GWWOBJ = ../pw4gww/fft_custom.o ../pw4gww/stop_pp.o ../pw4gww/mp_wave_parallel.o
GWWOBJ = ../gww/libgww.a ../minpack/minpacklib.a
@ -58,7 +56,7 @@ all : tldeps bse_main.x
bse_main.x : bse_main.o libbse.a $(BSEOBJS) $(PWOBJS) $(QEMODS) $(GWWOBJ)
$(LD) $(LDFLAGS) -o $@ \
bse_main.o libbse.a $(PWOBJS) $(QEMODS) $(LIBOBJS) $(QELIBS) $(PW4GWWOBJ) $(GWWOBJ)
bse_main.o libbse.a $(PWOBJS) $(QEMODS) $(QELIBS) $(PW4GWWOBJ) $(GWWOBJ)
- ( cd ../../bin ; ln -fs ../GWW/bse/$@ . )
tldeps:

12
GWW/gww/Makefile
View File

@ -43,9 +43,7 @@ times_gw.o \
vcprim.o
QEMODS = ../../Modules/libqemod.a ../../FFTXlib/libqefft.a \
../../KS_Solvers/libks_solvers.a ../../XClib/xc_lib.a \
../../LAXlib/libqela.a ../../UtilXlib/libutil.a ../../upflib/libupf.a
QEMODS = ../../KS_Solvers/libks_solvers.a $(BASEMODS)
LIBMIN= ../minpack/minpacklib.a
@ -53,14 +51,14 @@ TLDEPS=phlibs
all : tldeps gww.x gww_fit.x libgww.a
gww.x : gww.o $(GWWOBJS) $(LIBOBJS) $(QEMODS) $(LIBMIN)
gww.x : gww.o $(GWWOBJS) $(QEMODS) $(LIBMIN)
$(MPIF90) $(LDFLAGS) -o $@ \
gww.o $(GWWOBJS) $(QEMODS) $(LIBOBJS) $(LIBMIN) $(QELIBS)
gww.o $(GWWOBJS) $(QEMODS) $(LIBMIN) $(QELIBS)
- ( cd ../../bin; ln -fs ../GWW/gww/$@ . )
gww_fit.x : gww_fit.o $(GWWOBJS) $(LIBOBJS) $(QEMODS) $(LIBMIN)
gww_fit.x : gww_fit.o $(GWWOBJS) $(QEMODS) $(LIBMIN)
$(MPIF90) $(LDFLAGS) -o $@ \
gww_fit.o $(GWWOBJS) $(QEMODS) $(LIBOBJS) $(LIBMIN) $(QELIBS)
gww_fit.o $(GWWOBJS) $(QEMODS) $(LIBMIN) $(QELIBS)
- ( cd ../../bin; ln -fs ../GWW/gww/$@ . )
libgww.a : $(GWWOBJS)

13
GWW/head/Makefile
View File

@ -19,11 +19,10 @@ gwl_readin.o \
lanczos_k.o \
solve_head.o
QEMODS = ../../Modules/libqemod.a ../../upflib/libupf.a \
../../KS_Solvers/libks_solvers.a ../../XClib/xc_lib.a \
../../FFTXlib/libqefft.a ../../LAXlib/libqela.a ../../UtilXlib/libutil.a \
../../dft-d3/libdftd3qe.a
LIBPWPH = ../../PHonon/PH/libph.a ../..//LR_Modules/liblrmod.a ../../PW/src/libpw.a
QEMODS = $(BASEMODS)
LIBPWPH = ../../PHonon/PH/libph.a ../../LR_Modules/liblrmod.a \
../../PW/src/libpw.a ../../KS_Solvers/libks_solvers.a ../../dft-d3/libdftd3qe.a
LIBMIN= ../minpack/minpacklib.a
@ -31,9 +30,9 @@ TLDEPS= phlibs
all : tldeps head.x
head.x : head.o $(LIBOBJS) $(HEADOBJS) $(LIBPWPH) $(LIBMIN)
head.x : head.o $(HEADOBJS) $(LIBPWPH) $(LIBMIN)
$(LD) $(LDFLAGS) -o head.x head.o \
$(HEADOBJS) $(LIBPWPH) $(QEMODS) $(LIBOBJS) $(LIBMIN) $(QELIBS)
$(HEADOBJS) $(LIBPWPH) $(QEMODS) $(LIBMIN) $(QELIBS)
- ( cd ../../bin ; ln -fs ../GWW/head/head.x . )
tldeps :

4
GWW/head/close_gwl.f90
View File

@ -18,9 +18,9 @@ SUBROUTINE close_phq( flag )
!
USE io_global, ONLY : ionode, stdout
USE uspp, ONLY : okvan
USE units_ph, ONLY : iudwf, iubar, iudrhous, iuebar, iudrho, &
USE units_ph, ONLY : iubar, iudrhous, iuebar, iudrho, &
iudvscf, iucom, iudvkb3
USE units_lr, ONLY : iuwfc
USE units_lr, ONLY : iuwfc, iudwf
USE control_ph, ONLY : zue, epsil
USE recover_mod, ONLY : clean_recover
USE output, ONLY : fildrho, fildvscf

6
GWW/head/openfilq_gwl.f90
View File

@ -17,11 +17,11 @@ SUBROUTINE openfilq()
!
USE kinds, ONLY : DP
USE control_flags, ONLY : io_level, modenum
USE units_ph, ONLY : iudwf, iubar, iucom, iudvkb3, &
USE units_ph, ONLY : iubar, iucom, iudvkb3, &
iudrhous, iuebar, iudrho, iudyn, iudvscf, &
lrdwf, lrbar, lrcom, lrdvkb3, &
lrbar, lrcom, lrdvkb3, &
lrdrhous, lrebar, lrdrho
USE units_lr, ONLY : iuwfc, lrwfc
USE units_lr, ONLY : iuwfc, lrwfc, iudwf, lrdwf
! USE io_files, ONLY : tmp_dir
USE control_ph, ONLY : epsil, zue, ext_recover, trans, &
tmp_dir_ph, start_irr, last_irr

14
GWW/pw4gww/Makefile
View File

@ -49,23 +49,17 @@ gww_realus.o \
operator_1_vp.o \
operator_debug.o
QEMODS = ../../Modules/libqemod.a ../../upflib/libupf.a \
../../KS_Solvers/libks_solvers.a ../../XClib/xc_lib.a \
../../FFTXlib/libqefft.a ../../LAXlib/libqela.a \
../../UtilXlib/libutil.a ../../dft-d3/libdftd3qe.a
# dft-d3 required by xlf for obscure reasons
PWOBJS = ../../PW/src/libpw.a
QEMODS = $(BASEMODS)
PWOBJS = ../../PW/src/libpw.a ../../KS_Solvers/libks_solvers.a ../../dft-d3/libdftd3qe.a
PHOBJS = ../../PHonon/Gamma/libphcg.a
LIBMIN=
TLDEPS= phlibs
all : tldeps pw4gww.x
pw4gww.x : pw4gww.o libpw4gww.a $(PW4GWWOBJS) $(PWOBJS) $(QEMODS) $(LIBOBJS) $(LIBMIN)
pw4gww.x : pw4gww.o libpw4gww.a $(PW4GWWOBJS) $(PWOBJS) $(QEMODS)
$(LD) $(LDFLAGS) -o $@ \
pw4gww.o libpw4gww.a $(PWOBJS) $(PHOBJS) $(QEMODS) $(LIBOBJS) $(LIBMIN) $(QELIBS)
pw4gww.o libpw4gww.a $(PWOBJS) $(PHOBJS) $(QEMODS) $(QELIBS)
- ( cd ../../bin ; ln -fs ../GWW/pw4gww/$@ . )
tldeps :

11
GWW/simple/Makefile
View File

@ -22,13 +22,8 @@ SIMPLEOBJS = \
commutator.o
QEMODS = ../../Modules/libqemod.a ../../FFTXlib/libqefft.a \
../../upflib/libupf.a ../../KS_Solvers/libks_solvers.a \
../../LAXlib/libqela.a ../../UtilXlib/libutil.a \
../../dft-d3/libdftd3qe.a ../../XClib/xc_lib.a
# dft-d3 required by xlf for obscure reasons
PWOBJS = ../../PW/src/libpw.a
QEMODS = $(BASEMODS)
PWOBJS = ../../PW/src/libpw.a ../../KS_Solvers/libks_solvers.a ../../dft-d3/libdftd3qe.a
GWWOBJ = ../gww/libgww.a
TLDEPS= pwlibs gwwlib
@ -37,7 +32,7 @@ all : tldeps simple.x
simple.x : simple.o libsimple.a $(SIMPLEOBJS) $(PWOBJS) $(QEMODS) $(GWWOBJ)
$(LD) $(LDFLAGS) -o $@ \
simple.o libsimple.a $(PWOBJS) $(GWWOBJ) $(QEMODS) $(LIBOBJS) $(QELIBS) $(LIBMIN)
simple.o libsimple.a $(PWOBJS) $(GWWOBJ) $(QEMODS) $(QELIBS) $(LIBMIN)
- ( cd ../../bin ; ln -fs ../GWW/simple/$@ . )
tldeps :

8
GWW/simple_bse/Makefile
View File

@ -19,17 +19,15 @@ spectrum.o \
lanczos.o \
build_eemat.o
QEMODS = ../../Modules/libqemod.a ../../FFTXlib/libqefft.a \
../../KS_Solvers/libks_solvers.a \
../../LAXlib/libqela.a ../../UtilXlib/libutil.a ../../upflib/libupf.a
QEMODS = $(BASEMODS)
TLDEPS=bindir libs mods gwwlib
all : tldeps simple_bse.x libsimple_exc.a
simple_bse.x : simple_bse.o $(EXCOBJS) $(LIBOBJS) $(QEMODS)
simple_bse.x : simple_bse.o $(EXCOBJS) $(QEMODS)
$(MPIF90) $(LDFLAGS) -o $@ \
simple_bse.o $(EXCOBJS) $(QEMODS) $(LIBOBJS) $(QELIBS)
simple_bse.o $(EXCOBJS) $(QEMODS) $(QELIBS)
- ( cd ../../bin; ln -fs ../GWW/simple_bse/$@ . )

10
GWW/simple_ip/Makefile
View File

@ -16,10 +16,8 @@ diagonalization.o \
dielectric.o
QEMODS = ../../Modules/libqemod.a ../../FFTXlib/libqefft.a \
../../KS_Solvers/libks_solvers.a \
../../LAXlib/libqela.a ../../UtilXlib/libutil.a ../../upflib/libupf.a
PWOBJS = ../../PW/src/libpw.a
QEMODS = $(BASEMODS)
PWOBJS = ../../PW/src/libpw.a ../../KS_Solvers/libks_solvers.a
LIBMIN=
@ -27,9 +25,9 @@ TLDEPS=pwlibs
all : tldeps simple_ip.x libsimple_ip.a
simple_ip.x : simple_ip.o $(SIMPLEIPOBJS) $(LIBOBJS) $(PWOBJS) $(QEMODS) $(LIBMIN)
simple_ip.x : simple_ip.o $(SIMPLEIPOBJS) $(PWOBJS) $(QEMODS) $(LIBMIN)
$(MPIF90) $(LDFLAGS) -o $@ \
simple_ip.o $(SIMPLEIPOBJS) $(PWOBJS) $(QEMODS) $(LIBOBJS) $(LIBMIN) $(QELIBS)
simple_ip.o $(SIMPLEIPOBJS) $(PWOBJS) $(QEMODS) $(LIBMIN) $(QELIBS)
- ( cd ../../bin; ln -fs ../GWW/simple_ip/$@ . )
libsimple_ip.a : $(SIMPLEIPOBJS)

10
HP/src/Makefile
View File

@ -53,10 +53,8 @@ hp_write_chi_full.o \
hp_write_dnsq.o
LRMODS = ../../LR_Modules/liblrmod.a
PWOBJS = ../../PW/src/libpw.a
QEMODS = ../../Modules/libqemod.a ../../KS_Solvers/libks_solvers.a \
../../upflib/libupf.a ../../FFTXlib/libqefft.a ../../LAXlib/libqela.a\
../../UtilXlib/libutil.a ../../dft-d3/libdftd3qe.a ../../XClib/xc_lib.a
PWOBJS = ../../PW/src/libpw.a ../../KS_Solvers/libks_solvers.a ../../dft-d3/libdftd3qe.a
QEMODS = $(BASEMODS)
TLDEPS= hplibs
@ -64,9 +62,9 @@ all : tldeps libs-hp hp.x
libs-hp : libhp.a
hp.x : hp_main.o libhp.a $(PWOBJS) $(LRMODS) $(LIBOBJS)
hp.x : hp_main.o libhp.a $(PWOBJS) $(LRMODS)
$(LD) $(LDFLAGS) -o $@ hp_main.o libhp.a \
$(LRMODS) $(PWOBJS) $(QEMODS) $(LIBOBJS) $(QELIBS)
$(LRMODS) $(PWOBJS) $(QEMODS) $(QELIBS)
- ( cd ../../bin ; ln -fs ../HP/src/hp.x . )
tldeps :

6
HP/src/hp_allocate_q.f90
View File

@ -26,8 +26,7 @@ subroutine hp_allocate_q
USE eqv, ONLY : dpsi, evq, dmuxc, dvpsi
USE control_lr, ONLY : lgamma
USE ldaU, ONLY : Hubbard_lmax, nwfcU
USE ldaU_hp, ONLY : this_pert_is_on_file, &
swfcatomk, swfcatomkpq
USE ldaU_lr, ONLY : swfcatomk, swfcatomkpq
!
IMPLICIT NONE
INTEGER :: ik
@ -44,9 +43,6 @@ subroutine hp_allocate_q
ALLOCATE (dpsi(npwx*npol,nbnd))
ALLOCATE (dmuxc(dfftp%nnr,nspin_mag,nspin_mag))
!
ALLOCATE (this_pert_is_on_file(nksq))
this_pert_is_on_file(:) = .FALSE.
!
IF (okvan) THEN
ALLOCATE (eigqts(nat))
ALLOCATE (becp1(nksq))

8
HP/src/hp_close_q.f90
View File

@ -16,8 +16,8 @@ SUBROUTINE hp_close_q ( flag )
!
USE buffers, ONLY : close_buffer
USE io_files, ONLY : iunhub
USE units_lr, ONLY : iuwfc, iuatswfc
USE ldaU_hp, ONLY : iudwfc, iudvwfc
USE units_lr, ONLY : iuwfc, iuatswfc, iudwf
USE ldaU_hp, ONLY : iudvwfc
USE control_lr, ONLY : lgamma
!
IMPLICIT NONE
@ -28,10 +28,10 @@ SUBROUTINE hp_close_q ( flag )
CALL close_buffer(iuwfc,'delete')
!
IF (flag) THEN
CALL close_buffer(iudwfc,'delete')
CALL close_buffer(iudwf,'delete')
CALL close_buffer(iudvwfc,'delete')
ELSE
CALL close_buffer(iudwfc,'keep')
CALL close_buffer(iudwf,'keep')
CALL close_buffer(iudvwfc,'keep')
ENDIF
!

3
HP/src/hp_dealloc_1.f90
View File

@ -14,7 +14,8 @@ SUBROUTINE hp_dealloc_1()
!
USE lr_symm_base, ONLY : rtau
USE start_k, ONLY : xk_start, wk_start
USE ldaU_hp, ONLY : Rvect, dnsscf, dns0, dnsscf_tot, dns0_tot, &
USE ldaU_lr, ONLY : dnsscf
USE ldaU_hp, ONLY : Rvect, dns0, dnsscf_tot, dns0_tot, &
x_q, comp_iq
!
IMPLICIT NONE

6
HP/src/hp_dealloc_q.f90
View File

@ -22,8 +22,7 @@ SUBROUTINE hp_dealloc_q()
& dvxc_s, vsgga, segni
USE eqv, ONLY : dmuxc, dpsi, dvpsi, evq
USE control_lr, ONLY : lgamma, nbnd_occ
USE ldaU_hp, ONLY : this_pert_is_on_file, &
swfcatomk, swfcatomkpq
USE ldaU_lr, ONLY : swfcatomk, swfcatomkpq
!
IMPLICIT NONE
INTEGER :: ik
@ -42,9 +41,6 @@ SUBROUTINE hp_dealloc_q()
if (allocated(ikqs)) deallocate (ikqs)
if (allocated(m_loc)) deallocate (m_loc)
!
if (allocated(this_pert_is_on_file)) &
& deallocate (this_pert_is_on_file)
!
IF (okvan) THEN
if (allocated(eigqts)) deallocate (eigqts)
if (allocated(becp1)) then

8
HP/src/hp_dnsq.f90
View File

@ -38,12 +38,12 @@ SUBROUTINE hp_dnsq (lmetq0, iter, conv_root, dnsq)
USE control_flags, ONLY : iverbosity
USE qpoint, ONLY : nksq, ikks, ikqs
USE control_lr, ONLY : lgamma, nbnd_occ
USE units_lr, ONLY : iuwfc, lrwfc, iuatswfc
USE units_lr, ONLY : iuwfc, lrwfc, iuatswfc, iudwf, lrdwf
USE lr_symm_base, ONLY : nsymq
USE ldaU, ONLY : Hubbard_lmax, Hubbard_l, is_hubbard, offsetU, nwfcU
USE ldaU_hp, ONLY : conv_thr_chi, trace_dns_tot_old, &
conv_thr_chi_best, iter_best, iudwfc, lrdwfc, &
swfcatomk, swfcatomkpq
conv_thr_chi_best, iter_best
USE ldaU_lr, ONLY : swfcatomk, swfcatomkpq
USE hp_efermi_shift, ONLY : def
!
IMPLICIT NONE
@ -120,7 +120,7 @@ SUBROUTINE hp_dnsq (lmetq0, iter, conv_root, dnsq)
!
! At each SCF iteration for each ik read dpsi from file
!
CALL get_buffer (dpsi, lrdwfc, iudwfc, ik)
CALL get_buffer (dpsi, lrdwf, iudwf, ik)
!
! Loop on Hubbard atoms
!

3
HP/src/hp_dnstot_sum_q.f90
View File

@ -29,7 +29,8 @@ SUBROUTINE hp_dnstot_sum_q
USE control_flags, ONLY : iverbosity
USE lr_symm_base, ONLY : nsymq, invsymq, minus_q, rtau
USE ldaU, ONLY : Hubbard_lmax, Hubbard_l, is_hubbard
USE ldaU_hp, ONLY : nqsh, Rvect, dnsscf, dns0, dnsscf_tot, dns0_tot, &
USE ldaU_lr, ONLY : dnsscf
USE ldaU_hp, ONLY : nqsh, Rvect, dns0, dnsscf_tot, dns0_tot, &
skip_equivalence_q, nq1, nq2, nq3, x_q, nqs
!
IMPLICIT NONE

23
HP/src/hp_dvpsi_pert.f90
View File

@ -20,8 +20,6 @@ subroutine hp_dvpsi_pert (ik)
!
! dvpsi is for a given "k", "q" and "J"
!
! dvpsi is READ from file if this_pert_is_on_file(ik) = .TRUE.
! otherwise dvpsi is COMPUTED and WRITTEN on file
! (evc, swfcatomk, swfcatomkpq must be set)
!
USE kinds, ONLY : DP
@ -38,8 +36,8 @@ subroutine hp_dvpsi_pert (ik)
USE units_lr, ONLY : iuatswfc
USE control_lr, ONLY : lgamma
USE ldaU, ONLY : Hubbard_lmax, Hubbard_l, offsetU, nwfcU
USE ldaU_hp, ONLY : nqsh, perturbed_atom, this_pert_is_on_file, &
iudvwfc, lrdvwfc, swfcatomk, swfcatomkpq
USE ldaU_hp, ONLY : nqsh, perturbed_atom, iudvwfc, lrdvwfc
USE ldaU_lr, ONLY : swfcatomk, swfcatomkpq
!
IMPLICIT NONE
!
@ -59,23 +57,11 @@ subroutine hp_dvpsi_pert (ik)
DO na = 1, nat
IF (perturbed_atom(na)) counter = counter + 1
ENDDO
IF (counter.NE.1) CALL errore( 'hp_dvpsi_pert', "One perturbed atom must be specified", 1)
IF (counter /= 1) CALL errore( 'hp_dvpsi_pert', "One perturbed atom must be specified", 1)
!
dvpsi(:,:) = (0.0d0, 0.0d0)
!
! If this is not the first iteration, hence dvpsi was already
! computed before. So read it from file and exit.
!
IF (this_pert_is_on_file(ik)) THEN
!
CALL get_buffer(dvpsi, lrdvwfc, iudvwfc, ik)
CALL stop_clock ('hp_dvpsi_pert')
RETURN
!
ENDIF
!
! If this is a first iteration, then dvpsi must be computed
! and written on file.
! Compute dvpsi for ik and write on buffer iudvwfc
!
ALLOCATE (proj(nbnd,nwfcU))
!
@ -126,7 +112,6 @@ subroutine hp_dvpsi_pert (ik)
! Write dvpsi on file.
!
CALL save_buffer(dvpsi, lrdvwfc, iudvwfc, ik)
this_pert_is_on_file(ik) = .true.
!
DEALLOCATE (proj)
!

3
HP/src/hp_generate_grids.f90
View File

@ -20,7 +20,8 @@ SUBROUTINE hp_generate_grids()
USE qpoint, ONLY : xq
USE lr_symm_base, ONLY : rtau
USE ldaU, ONLY : Hubbard_lmax
USE ldaU_hp, ONLY : dnsscf, dns0, dnsscf_tot, dns0_tot, nqsh, &
USE ldaU_lr, ONLY : dnsscf
USE ldaU_hp, ONLY : dns0, dnsscf_tot, dns0_tot, nqsh, &
start_q, last_q, tmp_dir_hp, nqs, comp_iq
!
IMPLICIT NONE

10
HP/src/hp_openfil_q.f90
View File

@ -19,10 +19,10 @@ SUBROUTINE hp_openfil_q()
USE buffers, ONLY : open_buffer
USE qpoint, ONLY : nksq
USE control_lr, ONLY : lgamma
USE units_lr, ONLY : iuwfc, lrwfc, iuatswfc
USE units_lr, ONLY : iuwfc, lrwfc, iuatswfc, iudwf, lrdwf
USE ldaU, ONLY : nwfcU
USE ldaU_hp, ONLY : recalc_sym, tmp_dir_save, tmp_dir_hp, &
iudwfc, lrdwfc, iudvwfc, lrdvwfc
iudvwfc, lrdvwfc
!
IMPLICIT NONE
LOGICAL :: exst, exst_mem
@ -58,9 +58,9 @@ SUBROUTINE hp_openfil_q()
!
! Open a file to write/read a solution of the linear system (dpsi)
!
iudwfc = 22
lrdwfc = nbnd * npwx * npol
CALL open_buffer (iudwfc, 'dwfc', lrdwfc, io_level, exst_mem, exst, tmp_dir)
iudwf = 22
lrdwf = nbnd * npwx * npol
CALL open_buffer (iudwf, 'dwfc', lrdwf, io_level, exst_mem, exst, tmp_dir)
!
! Open a file to write/read S*phi at k and k+q (atomic wfct's)
!

3
HP/src/hp_read_dnsq.f90
View File

@ -18,7 +18,8 @@ SUBROUTINE hp_read_dnsq()
USE io_files, ONLY : prefix, tmp_dir
USE lsda_mod, ONLY : nspin
USE ldaU, ONLY : Hubbard_lmax, Hubbard_l, is_hubbard
USE ldaU_hp, ONLY : nah_pert, dns0, dnsscf, nqs, tmp_dir_hp
USE ldaU_lr, ONLY : dnsscf
USE ldaU_hp, ONLY : nah_pert, dns0, nqs, tmp_dir_hp
!
IMPLICIT NONE
!

273
HP/src/hp_solve_linear_system.f90
View File

@ -20,15 +20,13 @@ SUBROUTINE hp_solve_linear_system (na, iq)
USE io_global, ONLY : stdout
USE check_stop, ONLY : check_stop_now
USE wavefunctions, ONLY : evc
USE cell_base, ONLY : tpiba2
USE klist, ONLY : lgauss, ltetra, xk, wk, nelec, ngk, igk_k
USE gvect, ONLY : g
USE klist, ONLY : lgauss, ltetra, nelec, ngk
USE gvecs, ONLY : doublegrid
USE scf, ONLY : rho
USE fft_base, ONLY : dfftp, dffts
USE lsda_mod, ONLY : lsda, current_spin, isk
USE wvfct, ONLY : nbnd, npwx, g2kin, et
USE uspp, ONLY : okvan, vkb, nkb
USE wvfct, ONLY : nbnd, npwx
USE uspp, ONLY : okvan, nkb
USE uspp_param, ONLY : nhm
USE becmod, ONLY : allocate_bec_type, deallocate_bec_type, becp
USE buffers, ONLY : save_buffer, get_buffer
@ -39,20 +37,20 @@ SUBROUTINE hp_solve_linear_system (na, iq)
USE mp_pools, ONLY : inter_pool_comm, intra_pool_comm
USE mp, ONLY : mp_sum
USE hp_efermi_shift, ONLY : hp_ef_shift, def
USE eqv, ONLY : dvpsi, dpsi, evq
USE qpoint, ONLY : nksq, ikks, ikqs, xq
USE control_lr, ONLY : lgamma, nbnd_occ
USE qpoint, ONLY : nksq, ikks, xq
USE control_lr, ONLY : lgamma
USE units_lr, ONLY : iuwfc, lrwfc
USE lrus, ONLY : int3, int3_paw
USE dv_of_drho_lr, ONLY : dv_of_drho
USE fft_helper_subroutines
USE fft_interfaces, ONLY : fft_interpolate
USE lr_symm_base, ONLY : irotmq, minus_q, nsymq, rtau
USE ldaU_hp, ONLY : thresh_init, dnsscf, dns0, trace_dns_tot_old, &
USE ldaU_lr, ONLY : dnsscf
USE ldaU_hp, ONLY : thresh_init, dns0, trace_dns_tot_old, &
conv_thr_chi_best, iter_best, niter_max, nmix, &
alpha_mix, iudwfc, lrdwfc, code
USE apply_dpot_mod, ONLY : apply_dpot_allocate, apply_dpot_deallocate, &
apply_dpot_bands
alpha_mix, code, lrdvwfc, iudvwfc
USE apply_dpot_mod, ONLY : apply_dpot_allocate, apply_dpot_deallocate
USE response_kernels, ONLY : sternheimer_kernel
USE uspp_init, ONLY : init_us_2
!
IMPLICIT NONE
@ -63,20 +61,18 @@ SUBROUTINE hp_solve_linear_system (na, iq)
REAL(DP), ALLOCATABLE :: h_diag (:,:) ! diagonal part of the Hamiltonian
!
REAL(DP) :: thresh, & ! convergence threshold
anorm, & ! the norm of the error
averlt, & ! average number of iterations
dr2 ! self-consistency error
!
REAL(DP) :: dos_ef, & ! density of states at the Fermi level
weight, & ! Misc variables for metals
aux_avg(2) ! Misc variables for metals
REAL(DP) :: dos_ef
!! density of states at the Fermi level
!
REAL(DP), ALLOCATABLE :: becsum1(:,:,:)
!
COMPLEX(DP), ALLOCATABLE, TARGET :: dvscfin(:,:)
! change of the scf potential (input)
!
COMPLEX(DP), POINTER :: dvscfins(:,:)
COMPLEX(DP), POINTER :: dvscfins(:,:,:)
! change of the scf potential (smooth part only)
!
COMPLEX(DP), ALLOCATABLE :: drhoscf (:,:), &
@ -89,42 +85,27 @@ SUBROUTINE hp_solve_linear_system (na, iq)
ldoss (:,:), & ! as above, without augmentation charges
dbecsum (:,:,:,:), & ! the derivative of becsum
aux2 (:,:), & ! auxiliary arrays
mixin(:), mixout(:), & ! auxiliary arrays for mixing of the response potential
tg_dv(:,:), & ! Task groups: auxiliary array for potential * wfct
tg_psic(:,:) ! Task groups: auxiliary array for wavefunctions
mixin(:), mixout(:) ! auxiliary arrays for mixing of the response potential
COMPLEX(DP), ALLOCATABLE :: t(:,:,:,:), tmq(:,:,:)
! PAW: auxiliary arrays
LOGICAL :: conv_root, & ! true if linear system is converged
exst, & ! used to open the recover file
lmetq0, & ! true if xq=(0,0,0) in a metal
LOGICAL :: all_conv
!! True if sternheimer_kernel is converged at all k points
LOGICAL :: lmetq0, & ! true if xq=(0,0,0) in a metal
convt, & ! not needed for HP
convt_chi ! used instead of convt to control the convergence
REAL(DP), PARAMETER :: tr2 = 1.D-30 ! threshold parameter
INTEGER :: ibnd, & ! counter on bands
iter, & ! counter on iterations
lter, & ! counter on iterations of linear system
ltaver, & ! average counter
lintercall, & ! average number of calls to cgsolve_all
INTEGER :: iter, & ! counter on iterations
ik, ikk, & ! counter on k points
ikq, & ! counter on k+q points
ig, & ! counter on G vectors
ndim, &
is, & ! counter on spin polarizations
nt, & ! counter on types
ios, & ! integer variable for I/O control
incr, & ! used for task groups
v_siz, & ! size of the potential
npw, & ! number of plane waves at k
npwq ! number of plane waves at k+q
npw ! number of plane waves at k
REAL(DP) :: tcpu, get_clock ! timing variables
CHARACTER(LEN=256) :: filename, &
flmixdpot = 'mixd'
EXTERNAL ch_psi_all, cg_psi
CHARACTER(LEN=256) :: flmixdpot = 'mixd'
!
CALL start_clock ('hp_solve_linear_system')
!
@ -134,15 +115,16 @@ SUBROUTINE hp_solve_linear_system (na, iq)
!
! Allocate arrays for the SCF density/potential
!
ALLOCATE (drhoscf (dffts%nnr, nspin_mag))
ALLOCATE (drhoscf (dfftp%nnr, nspin_mag))
ALLOCATE (drhoscfh(dfftp%nnr, nspin_mag))
ALLOCATE (dvscfin (dfftp%nnr, nspin_mag))
ALLOCATE (dvscfout(dfftp%nnr, nspin_mag))
!
dvscfin = (0.0_DP, 0.0_DP)
IF (doublegrid) THEN
ALLOCATE (dvscfins(dffts%nnr, nspin_mag))
ALLOCATE (dvscfins(dffts%nnr, nspin_mag, 1))
ELSE
dvscfins => dvscfin
dvscfins(1:dffts%nnr, 1:nspin_mag, 1:1) => dvscfin
ENDIF
!
! USPP-specific allocations
@ -176,16 +158,6 @@ SUBROUTINE hp_solve_linear_system (na, iq)
convt = .FALSE.
convt_chi = .FALSE.
!
incr = 1
IF ( dffts%has_task_groups ) THEN
!
v_siz = dffts%nnr_tg
ALLOCATE( tg_dv ( v_siz, nspin_mag ) )
ALLOCATE( tg_psic( v_siz, npol ) )
incr = fftx_ntgrp(dffts)
!
ENDIF
!
! If q=0 for a metal: allocate and compute local DOS and DOS at Ef
!
lmetq0 = (lgauss .OR. ltetra) .AND. lgamma
@ -198,158 +170,65 @@ SUBROUTINE hp_solve_linear_system (na, iq)
IF (.NOT.okpaw) DEALLOCATE (becsum1)
ENDIF
!
! Compute dV_bare * psi and write to buffer iubar
!
DO ik = 1, nksq
!
ikk = ikks(ik)
npw = ngk(ikk)
!
IF (lsda) current_spin = isk(ikk)
!
! Read unperturbed KS wavefuctions psi(k) and psi(k+q)
!
IF (nksq > 1) THEN
CALL get_buffer(evc, lrwfc, iuwfc, ikk)
ENDIF
!
! Computes (iter=1) or reads (iter>1) the action of the perturbing
! potential on the unperturbed KS wavefunctions: |dvpsi> = dV_pert * |evc>
! See Eq. (46) in Ref. [1]
!
CALL hp_dvpsi_pert(ik)
!
ENDDO ! ik
!
! The loop of the linear-response calculation
!
DO iter = 1, niter_max
!
WRITE(stdout,'(/6x,"atom #",i3,3x,"q point #",i4,3x,"iter # ",i3)') na, iq, iter
!
ltaver = 0
lintercall = 0
!
drhoscf(:,:) = (0.d0, 0.d0)
dvscfout(:,:) = (0.d0, 0.d0)
dbecsum(:,:,:,:) = (0.d0, 0.d0)
!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!!!!!!!!!!!!!!!!! START OF THE K LOOP !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
IF ( iter == 1 ) THEN
! Starting threshold for iterative solution of the linear system.
! A strickt threshold for the first iteration is needed,
! because we need dns0 to very high precision.
thresh = thresh_init * nelec
ELSE
! Threshold for iterative solution of the linear system.
! We start with not a strict threshold for iter=2, and then
! it decreases with iterations.
thresh = MIN (1.D-1 * SQRT(dr2), 1.D-2)
ENDIF
!
DO ik = 1, nksq
!
ikk = ikks(ik)
ikq = ikqs(ik)
npw = ngk(ikk)
npwq = ngk(ikq)
!
IF (lsda) current_spin = isk(ikk)
!
! Read unperturbed KS wavefuctions psi(k) and psi(k+q)
!
IF (nksq.gt.1) THEN
IF (lgamma) THEN
CALL get_buffer (evc, lrwfc, iuwfc, ikk)
ELSE
CALL get_buffer (evc, lrwfc, iuwfc, ikk)
CALL get_buffer (evq, lrwfc, iuwfc, ikq)
ENDIF
ENDIF
!
! USPP: Compute the projectors vkb at k+q
!
CALL init_us_2 (npwq, igk_k(1,ikq), xk(1,ikq), vkb)
!
! Compute the kinetic energy at k+q
!
CALL g2_kin (ikq)
!
! Compute preconditioning matrix h_diag used by cgsolve_all
!
CALL h_prec (ik, evq, h_diag)
!
! Computes (iter=1) or reads (iter>1) the action of the perturbing
! potential on the unperturbed KS wavefunctions: |dvpsi> = dV_pert * |evc>
! See Eq. (46) in Ref. [1]
!
CALL hp_dvpsi_pert(ik)
!
IF ( iter > 1 ) THEN
!
! Add the contribution of the self consistent term.
! Calculates dvscf_q*psi(k) in G-space, for all bands, k=ik
! dvscf_q from previous iteration (mix_potential)
!
CALL start_clock ('hp_vpsifft')
CALL apply_dpot_bands(ik, nbnd_occ(ikk), dvscfins, evc, aux2)
dvpsi = dvpsi + aux2
CALL stop_clock ('hp_vpsifft')
!
! USPP: there is an additional self-consistent term proportional to int3
! |dvpsi> = |dvpsi> + dV_HXC*|evc> + int3 * |beta><beta|evc>
!
IF (okvan) CALL adddvscf(1, ik)
!
ENDIF
!
! Ortogonalize dvpsi to valence states: ps = <evq|dvpsi>
! Apply -P_c^+. See Eq. (A21) in Ref. [1]
!
CALL orthogonalize(dvpsi, evq, ikk, ikq, dpsi, npwq, .FALSE.)
!
IF ( iter == 1 ) THEN
!
! At the first iteration dpsi and dvscfin are set to zero
!
dpsi(:,:) = (0.d0, 0.d0)
dvscfin(:,:) = (0.d0, 0.d0)
!
! Starting threshold for iterative solution of the linear system.
! A strickt threshold for the first iteration is needed,
! because we need dns0 to very high precision.
!
thresh = thresh_init * nelec
!
ELSE
!
! Starting value for dpsi is read from file
!
CALL get_buffer( dpsi, lrdwfc, iudwfc, ik)
!
! Threshold for iterative solution of the linear system.
! We start with not a strict threshold for iter=2, and then
! it decreases with iterations.
!
thresh = MIN (1.D-1 * SQRT(dr2), 1.D-2)
!
ENDIF
!
! Iterative solution of the linear system:
! (H + Q - eS) * |dpsi> = |dvpsi>,
! where |dvpsi> = - P_c^+ (dV_HXC + dV_pert) * |evc>
! See Eq. (43) in Ref. [1]
!
CALL cgsolve_all (ch_psi_all, cg_psi, et(1,ikk), dvpsi, dpsi, h_diag, &
& npwx, npwq, thresh, ik, lter, conv_root, anorm, nbnd_occ(ikk), npol )
!
ltaver = ltaver + lter
!
lintercall = lintercall + 1
!
IF (.NOT.conv_root) THEN
WRITE( stdout, '(6x,"kpoint",i4, &
& " hp_solve_linear_system: root not converged, thresh < ",e10.3)') ik , anorm
IF (iter == 1) WRITE( stdout, '(6x,"Try to increase thresh_init...")')
ENDIF
!
! Writes dpsi on file for a given k
!
CALL save_buffer (dpsi, lrdwfc, iudwfc, ik)
!
! Setup the weight at point k (normalized by the number of k points)
!
weight = wk(ikk)
!
! Calculates the response charge density (sum over k)
! See Eq. (48) in Ref. [1]
!
CALL incdrhoscf (drhoscf(:,current_spin), weight, ik, &
& dbecsum(:,:,current_spin,1), dpsi)
!
ENDDO ! k points
! Compute drhoscf, the charge density response to the total potential
!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!!!!!!!!!!!!!!!!!!! END OF THE K LOOP !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
CALL sternheimer_kernel(iter==1, .FALSE., 1, lrdvwfc, iudvwfc, &
thresh, dvscfins, all_conv, averlt, drhoscf, dbecsum, exclude_hubbard=.TRUE.)
!
#if defined (__MPI)
IF ((.NOT. all_conv) .AND. (iter == 1)) THEN
WRITE(stdout, '(6x, "sternheimer_kernel not converged. Try to increase thresh_init.")')
ENDIF
!
! USPP: The calculation of dbecsum is distributed across processors (see addusdbec)
! Sum over processors the contributions coming from each slice of bands
!
CALL mp_sum ( dbecsum, intra_pool_comm )
!
#endif
!
! Copy/interpolate the response density drhoscf -> drhoscfh
!
IF (doublegrid) THEN
@ -364,10 +243,9 @@ SUBROUTINE hp_solve_linear_system (na, iq)
!
IF (okvan) CALL lr_addusddens (drhoscfh, dbecsum)
!
#if defined (__MPI)
call mp_sum ( drhoscf, inter_pool_comm )
CALL mp_sum ( drhoscfh, inter_pool_comm )
IF (okpaw) CALL mp_sum ( dbecsum, inter_pool_comm )
#endif
!
! PAW: the factor of 2 is due to the presence of the CC term
! (see first two terms in Eq.(9) in PRB 81, 075123 (2010))
@ -452,7 +330,7 @@ SUBROUTINE hp_solve_linear_system (na, iq)
!
IF (doublegrid) THEN
DO is = 1, nspin_mag
CALL fft_interpolate (dfftp, dvscfin(:,is), dffts, dvscfins(:,is))
CALL fft_interpolate (dfftp, dvscfin(:,is), dffts, dvscfins(:,is,1))
ENDDO
ENDIF
!
@ -479,16 +357,7 @@ SUBROUTINE hp_solve_linear_system (na, iq)
!
IF ( iter == 1 ) dns0(:,:,:,:,iq) = dnsscf(:,:,:,:,iq)
!
! Compute the average number of iterations
!
#if defined (__MPI)
aux_avg(1) = DBLE(ltaver)
aux_avg(2) = DBLE(lintercall)
CALL mp_sum ( aux_avg, inter_pool_comm )
averlt = aux_avg(1) / aux_avg(2)
#else
averlt = DBLE(ltaver) / DBLE(lintercall)
#endif
! Print the average number of iterations
!
tcpu = get_clock(code)
!
@ -507,12 +376,10 @@ SUBROUTINE hp_solve_linear_system (na, iq)
CALL hp_stop_smoothly (.TRUE.)
ENDIF
!
IF (convt_chi) goto 155
IF (convt_chi) EXIT
!
ENDDO ! loop over the iterations iter
!
155 CONTINUE
!
CALL apply_dpot_deallocate()
DEALLOCATE (h_diag)
DEALLOCATE (aux2)
@ -534,10 +401,6 @@ SUBROUTINE hp_solve_linear_system (na, iq)
DEALLOCATE (tmq)
ENDIF
CALL deallocate_bec_type (becp)
IF ( dffts%has_task_groups ) THEN
DEALLOCATE( tg_dv )
DEALLOCATE( tg_psic )
ENDIF
!
WRITE( stdout,*) " "
WRITE( stdout,*) " =--------------------------------------------="

3
HP/src/hp_write_dnsq.f90
View File

@ -18,7 +18,8 @@ SUBROUTINE hp_write_dnsq(iq)
USE lsda_mod, ONLY : nspin
USE io_files, ONLY : prefix, tmp_dir
USE ldaU, ONLY : Hubbard_lmax, Hubbard_l, is_hubbard
USE ldaU_hp, ONLY : nah_pert, dns0, dnsscf, x_q
USE ldaU_lr, ONLY : dnsscf
USE ldaU_hp, ONLY : nah_pert, dns0, x_q
!
IMPLICIT NONE
!

7
HP/src/hpcom.f90
View File

@ -49,7 +49,6 @@ MODULE ldaU_hp
LOGICAL, ALLOCATABLE :: todo_atom(:), & ! Which atoms must be perturbed
perturbed_atom(:), & ! Controls which atom is perturbed in the HP
! calculation
this_pert_is_on_file(:), & ! The perturbation is written on file or not
comp_iq(:) ! If .true. this q point has to be calculated
!
INTEGER :: nath, & ! Number of (real) atoms in the primitive cell
@ -77,8 +76,6 @@ MODULE ldaU_hp
nq1, nq2, nq3, & ! Number of q points in each direction
nqs, & ! Number of q points to be calculated
start_q, last_q, & ! Initial and final q in the list
iudwfc, & ! Unit for response wavefunctions
lrdwfc, & ! Length of the record for response wavefunctions
iudvwfc, & ! Unit for the perturbing potential * wavefunctions
lrdvwfc ! Length of the record for the perturbing potential * wavefunctions
!
@ -116,7 +113,6 @@ MODULE ldaU_hp
!
COMPLEX(DP), ALLOCATABLE :: dns0(:,:,:,:,:), & ! Bare response occupation matrix
! (from 1st iteration)
dnsscf(:,:,:,:,:), & ! SCF response occupation matrix
dns0_tot(:,:,:,:,:), & ! Total bare response occupation matrix
! (summed over q)
dnsscf_tot(:,:,:,:,:), & ! Total SCF response occupation matrix
@ -126,7 +122,4 @@ MODULE ldaU_hp
!
INTEGER, ALLOCATABLE :: ityp_new(:) ! Types of atoms
!
COMPLEX(DP), ALLOCATABLE, TARGET :: swfcatomk(:,:) ! S * atomic wfc at k
COMPLEX(DP), POINTER :: swfcatomkpq(:,:) ! S * atomic wfc at k+q
!
END MODULE ldaU_hp

6
HP/src/make.depend
View File

@ -9,7 +9,6 @@ hp_allocate_q.o : ../../Modules/wavefunctions.o
hp_allocate_q.o : ../../PW/src/ldaU.o
hp_allocate_q.o : ../../PW/src/pwcom.o
hp_allocate_q.o : ../../upflib/uspp.o
hp_allocate_q.o : hpcom.o
hp_bcast_input.o : ../../LR_Modules/lrcom.o
hp_bcast_input.o : ../../Modules/check_stop.o
hp_bcast_input.o : ../../Modules/control_flags.o
@ -49,7 +48,6 @@ hp_dealloc_q.o : ../../LR_Modules/lrcom.o
hp_dealloc_q.o : ../../Modules/becmod.o
hp_dealloc_q.o : ../../Modules/noncol.o
hp_dealloc_q.o : ../../upflib/uspp.o
hp_dealloc_q.o : hpcom.o
hp_dnsq.o : ../../LR_Modules/lrcom.o
hp_dnsq.o : ../../Modules/constants.o
hp_dnsq.o : ../../Modules/control_flags.o
@ -196,6 +194,7 @@ hp_read_chi.o : ../../Modules/io_global.o
hp_read_chi.o : ../../Modules/ions_base.o
hp_read_chi.o : ../../Modules/kind.o
hp_read_chi.o : hpcom.o
hp_read_dnsq.o : ../../LR_Modules/lrcom.o
hp_read_dnsq.o : ../../Modules/io_files.o
hp_read_dnsq.o : ../../Modules/io_global.o
hp_read_dnsq.o : ../../Modules/ions_base.o
@ -261,8 +260,8 @@ hp_solve_linear_system.o : ../../FFTXlib/fft_interfaces.o
hp_solve_linear_system.o : ../../LR_Modules/apply_dpot_mod.o
hp_solve_linear_system.o : ../../LR_Modules/dv_of_drho.o
hp_solve_linear_system.o : ../../LR_Modules/lrcom.o
hp_solve_linear_system.o : ../../LR_Modules/response_kernels.o
hp_solve_linear_system.o : ../../Modules/becmod.o
hp_solve_linear_system.o : ../../Modules/cell_base.o
hp_solve_linear_system.o : ../../Modules/check_stop.o
hp_solve_linear_system.o : ../../Modules/fft_base.o
hp_solve_linear_system.o : ../../Modules/io_global.o
@ -335,6 +334,7 @@ hp_write_chi.o : hpcom.o
hp_write_chi_full.o : ../../Modules/io_files.o
hp_write_chi_full.o : ../../Modules/ions_base.o
hp_write_chi_full.o : hpcom.o
hp_write_dnsq.o : ../../LR_Modules/lrcom.o
hp_write_dnsq.o : ../../Modules/io_files.o
hp_write_dnsq.o : ../../Modules/io_global.o
hp_write_dnsq.o : ../../Modules/ions_base.o

2
LAXlib/tests/Makefile
View File

@ -23,7 +23,7 @@ common: $(UTILS)
%.x: %.o
$(LD) $(LDFLAGS) $< $(UTILS) -o $@ ../libqela.a \
../../UtilXlib/libutil.a $(LIBOBJS) $(QELIBS)
../../UtilXlib/libutil.a $(QELIBS)
clean :
- /bin/rm -f *.o *.a *.d *.i *~ *_tmp.f90 *.mod *.L *.x

5
LR_Modules/CMakeLists.txt
View File

@ -1,5 +1,7 @@
set(sources
adddvscf.f90
adddvscf_ph_mag.f90
adddvhubscf.f90
addusdbec.f90
addusdbec_nc.f90
apply_dpot.f90
@ -56,7 +58,8 @@ set(sources
ccg_psi.f90
ccgsolve_all.f90
ch_psi_all_complex.f90
addusddenseq.f90)
addusddenseq.f90
response_kernels.f90)
qe_add_library(qe_lr_modules ${sources})
target_link_libraries(qe_lr_modules

5
LR_Modules/Makefile
View File

@ -7,6 +7,8 @@ MODFLAGS=$(BASEMOD_FLAGS) $(MOD_FLAG)$(TOPDIR)/PW/src
LR_MODULES = \
adddvscf.o \
adddvscf_ph_mag.o \
adddvhubscf.o \
addusdbec.o \
addusdbec_nc.o \
apply_dpot.o \
@ -64,7 +66,8 @@ lr_dot_magnons.o \
lanczos_pseudohermitian.o \
lanczos_nonhermitian.o \
lanczos_pseudohermitian_c.o \
lanczos_nonhermitian_c.o
lanczos_nonhermitian_c.o \
response_kernels.o
TLDEPS=mods pwlibs

4
PHonon/PH/adddvhubscf.f90 → LR_Modules/adddvhubscf.f90
View File

@ -33,7 +33,7 @@ SUBROUTINE adddvhubscf (ipert, ik)
USE ions_base, ONLY : nat, ityp, ntyp => nsp
USE ldaU, ONLY : Hubbard_lmax, Hubbard_l, offsetU, is_hubbard, &
Hubbard_J0, nwfcU
USE ldaU_ph, ONLY : dnsscf, swfcatomk, swfcatomkpq, proj1, effU
USE ldaU_lr, ONLY : effU, swfcatomk, swfcatomkpq, dnsscf
USE wvfct, ONLY : npwx, nbnd
USE control_lr, ONLY : lgamma, nbnd_occ
USE lsda_mod, ONLY : lsda, nspin, current_spin, isk
@ -57,7 +57,7 @@ SUBROUTINE adddvhubscf (ipert, ik)
INTEGER :: npw, npwq, ikk, ikq, op_spin
INTEGER :: i, j, k, nt, l, ih, n, ig, ihubst, ihubst1, ihubst2, &
nah, m, m1, m2, ibnd, ldim
COMPLEX(DP), ALLOCATABLE :: dvhubscf(:,:), dvqi(:,:)
COMPLEX(DP), ALLOCATABLE :: dvhubscf(:,:), dvqi(:,:), proj1(:, :)
COMPLEX(DP), EXTERNAL :: ZDOTC
!
CALL start_clock ( 'adddvhubscf' )

1
PHonon/PH/adddvscf_ph_mag.f90 → LR_Modules/adddvscf_ph_mag.f90
View File

@ -25,7 +25,6 @@ subroutine adddvscf_ph_mag (ipert, ik, becp1)
USE klist, ONLY : ngk
USE noncollin_module, ONLY : noncolin, npol
USE becmod, ONLY : bec_type
! modules from phcom
USE lrus, ONLY : int3, int3_nc
USE qpoint, ONLY : nksq, ikks, ikqs
USE eqv, ONLY : dvpsi

23
LR_Modules/lrcom.f90
View File

@ -200,6 +200,27 @@ MODULE units_lr
INTEGER :: iuwfc, & ! unit for wavefunctions
lrwfc, & ! the length of wavefunction record
iuatwfc, & ! unit for atomic wavefunctions
iuatswfc ! unit for atomic wavefunctions * S
iuatswfc,& ! unit for atomic wavefunctions * S
iudwf, & ! unit with D psi
lrdwf ! length of D psi record
!
END MODULE units_lr
MODULE ldaU_lr
!
USE kinds, ONLY : DP
USE parameters, ONLY : ntypx
!
REAL(DP) :: effU(ntypx)
! effective Hubbard parameter: effU = Hubbard_U - Hubbard_J0
! TODO: Can be moved to PW/ldaU
!
COMPLEX(DP), ALLOCATABLE :: dnsscf(:,:,:,:,:)
!! SCF derivative of ns
!
COMPLEX(DP), ALLOCATABLE, TARGET :: swfcatomk(:,:)
!! S * atomic wfc at k
COMPLEX(DP), POINTER :: swfcatomkpq(:,:)
!! S * atomic wfc at k+q
!
END MODULE ldaU_lr

32
LR_Modules/make.depend
View File

@ -17,12 +17,29 @@ adddvepsi_us.o : ../Modules/noncol.o
adddvepsi_us.o : ../PW/src/pwcom.o
adddvepsi_us.o : ../upflib/uspp.o
adddvepsi_us.o : lrcom.o
adddvhubscf.o : ../Modules/io_files.o
adddvhubscf.o : ../Modules/ions_base.o
adddvhubscf.o : ../Modules/kind.o
adddvhubscf.o : ../Modules/mp_bands.o
adddvhubscf.o : ../Modules/wavefunctions.o
adddvhubscf.o : ../PW/src/buffers.o
adddvhubscf.o : ../PW/src/ldaU.o
adddvhubscf.o : ../PW/src/pwcom.o
adddvhubscf.o : ../UtilXlib/mp.o
adddvhubscf.o : lrcom.o
adddvscf.o : ../Modules/ions_base.o
adddvscf.o : ../Modules/kind.o
adddvscf.o : ../Modules/noncol.o
adddvscf.o : ../PW/src/pwcom.o
adddvscf.o : ../upflib/uspp.o
adddvscf.o : lrcom.o
adddvscf_ph_mag.o : ../Modules/becmod.o
adddvscf_ph_mag.o : ../Modules/ions_base.o
adddvscf_ph_mag.o : ../Modules/kind.o
adddvscf_ph_mag.o : ../Modules/noncol.o
adddvscf_ph_mag.o : ../PW/src/pwcom.o
adddvscf_ph_mag.o : ../upflib/uspp.o
adddvscf_ph_mag.o : lrcom.o
addusdbec.o : ../Modules/becmod.o
addusdbec.o : ../Modules/ions_base.o
addusdbec.o : ../Modules/kind.o
@ -321,6 +338,7 @@ lr_sym_mod.o : ../PW/src/symm_base.o
lr_sym_mod.o : lrcom.o
lrcom.o : ../Modules/becmod.o
lrcom.o : ../Modules/kind.o
lrcom.o : ../Modules/parameters.o
mix_pot.o : ../Modules/io_files.o
mix_pot.o : ../Modules/kind.o
mix_pot.o : ../Modules/mp_bands.o
@ -357,6 +375,20 @@ qdipol_cryst.o : ../Modules/kind.o
qdipol_cryst.o : ../PW/src/pwcom.o
qdipol_cryst.o : ../upflib/uspp.o
qdipol_cryst.o : lrcom.o
response_kernels.o : ../Modules/fft_base.o
response_kernels.o : ../Modules/io_global.o
response_kernels.o : ../Modules/ions_base.o
response_kernels.o : ../Modules/kind.o
response_kernels.o : ../Modules/mp_pools.o
response_kernels.o : ../Modules/noncol.o
response_kernels.o : ../Modules/wavefunctions.o
response_kernels.o : ../PW/src/buffers.o
response_kernels.o : ../PW/src/ldaU.o
response_kernels.o : ../PW/src/pwcom.o
response_kernels.o : ../UtilXlib/mp.o
response_kernels.o : ../upflib/uspp.o
response_kernels.o : apply_dpot_mod.o
response_kernels.o : lrcom.o
set_dbecsum_nc.o : ../Modules/ions_base.o
set_dbecsum_nc.o : ../Modules/kind.o
set_dbecsum_nc.o : ../Modules/noncol.o

284
LR_Modules/response_kernels.f90 Normal file
View File

@ -0,0 +1,284 @@
!
! Copyright (C) 2001-2018 Quantum ESPRESSO group
! This file is distributed under the terms of the
! GNU General Public License. See the file `License'
! in the root directory of the present distribution,
! or http://www.gnu.org/copyleft/gpl.txt .
!
!
!------------------------------------------------------------------------------
MODULE response_kernels
CONTAINS
SUBROUTINE sternheimer_kernel(first_iter, time_reversed, npert, lrdvpsi, iudvpsi, &
thresh, dvscfins, all_conv, avg_iter, drhoout, dbecsum, dbecsum_nc, &
exclude_hubbard)
!----------------------------------------------------------------------------
!! Compute the density response to the perturbation dV = dV_bare + dV_ind by the
!! non-interacting susceptibility. Solve Sternheimer equation
!! (H - e) * dpsi = dvpsi = -P_c+ * (dV_bare + dV_ind) * psi.
!!
!! Perfoms the following tasks:
!! a) reads the bare potential term Delta V | psi > from buffer iudvpsi.
!! b) adds to it the screening term Delta V_{SCF} | psi >.
!! If lda_plus_u=.true. compute also the SCF part
!! of the response Hubbard potential.
!! c) applies P_c^+ (orthogonalization to valence states).
!! d) calls cgsolve_all to solve the linear system.
!! e) returns the Delta rho, and if lda_plus_u=.true. also return dbecsum
!!
!! dV_bare * psi is read from buffer iudvpsi, so they must be already calculated.
!! dV_ind is given by input dvscfins, and dV_ind * psi is calculated in apply_dpot_bands.
!!
!! For USPPs, adddvscf is called, so relevant arrays must be already calculated.
!! For DFT+U, adddvhubscf is called, so relevant arrays must be already calculated.
!!
!! Results are added to drhoout, dbecsum, dbecsum_nc, so these output arrays should
!! be initialized before calling this subroutine.
!!
!! Input:
!! - first_iter: true if the first iteration, where dvscfins = 0
!! - time_reversed: false for ordinary cases, true for time-reversed wave functions
!! which is need in the noncolliner magnetic case
!! - npert: number of perturbations
!! - lrdvpsi: record length for the buffer storing dV_bare * psi
!! - iudvpsi: unit for the buffer storing dV_bare * psi
!! - thresh: threshold for solving Sternheimer equation
!! - dvscfins: dV_ind calculated in the previous iteration
!! - exclude_hubbard: If TRUE, do not add the response of the Hubbard potential.
!! Used in hp.x (Default: FALSE)
!!
!! Output:
!! - avg_iter: average number of iterations for the linear equation solver
!! - drhoout: induced charge density
!! - dbecsum: becsum with dpsi
!! - dbecsum_nc: becsum with dpsi. Optional, used if noncolin is true.
!!
!! FIXME: Only 1:dffts%nnr is used for drhoout, but it is allocated as 1:dfftp%nnr.
!----------------------------------------------------------------------------
USE kinds, ONLY : DP
USE io_global, ONLY : stdout
USE mp, ONLY : mp_sum
USE mp_pools, ONLY : inter_pool_comm
USE buffers, ONLY : get_buffer, save_buffer
USE fft_base, ONLY : dfftp
USE ions_base, ONLY : nat
USE klist, ONLY : xk, wk, ngk, igk_k
USE lsda_mod, ONLY : lsda, nspin, current_spin, isk
USE spin_orb, ONLY : domag
USE wvfct, ONLY : nbnd, npwx, et
USE wavefunctions, ONLY : evc
USE noncollin_module, ONLY : noncolin, npol, nspin_mag
USE uspp, ONLY : vkb
USE uspp_param, ONLY : nhm
USE ldaU, ONLY : lda_plus_u
USE units_lr, ONLY : iuwfc, lrwfc, lrdwf, iudwf
USE control_lr, ONLY : nbnd_occ, lgamma
USE qpoint, ONLY : nksq, ikks, ikqs
USE qpoint_aux, ONLY : ikmks, ikmkmqs, becpt
USE eqv, ONLY : dpsi, dvpsi, evq
USE apply_dpot_mod, ONLY : apply_dpot_bands
!
IMPLICIT NONE
!
LOGICAL, INTENT(IN) :: first_iter
!! true if the first iteration.
LOGICAL, INTENT(IN) :: time_reversed
!! true if solving for time reversed wave functions
LOGICAL, INTENT(IN), OPTIONAL :: exclude_hubbard
!! true if ignoring the Hubbard response term
INTEGER, INTENT(IN) :: npert
!! number of perturbations
INTEGER, INTENT(IN) :: lrdvpsi
!! record length for the buffer storing dV_bare * psi
INTEGER, INTENT(IN) :: iudvpsi
!! unit for the buffer storing dV_bare * psi
REAL(DP), INTENT(IN) :: thresh
!! threshold for solving the linear equation
LOGICAL, INTENT(OUT) :: all_conv
!! True if converged at all k points and perturbations
REAL(DP), INTENT(OUT) :: avg_iter
!! average number of iterations for the linear equation solver
COMPLEX(DP), POINTER, INTENT(IN) :: dvscfins(:, :, :)
!! dV_ind calculated in the previous iteration
COMPLEX(DP), INTENT(INOUT) :: drhoout(dfftp%nnr, nspin_mag, npert)
!! induced charge density
COMPLEX(DP), INTENT(INOUT) :: dbecsum(nhm*(nhm+1)/2, nat, nspin_mag, npert)
!! becsum with dpsi
COMPLEX(DP), INTENT(INOUT), OPTIONAL :: dbecsum_nc(nhm, nhm, nat, nspin, npert)
!! becsum with dpsi. Used if noncolin is true.
!
LOGICAL :: conv_root
!! true if linear system is converged
LOGICAL :: exclude_hubbard_
!! Local variable to set the default of exclude_hubbard to false
INTEGER :: ikk, ikq, npw, npwq, ipert, num_iter, ik, nrec, ikmk, ikmkmq
!! counters
INTEGER :: tot_num_iter
!! total number of iterations in cgsolve_all
INTEGER :: tot_cg_calls
!! total number of cgsolve_all calls
REAL(DP) :: anorm
!! the norm of the error of the conjugate gradient solution
REAL(DP) :: rsign
!! sign of the term in the magnetization
REAL(DP), ALLOCATABLE :: h_diag(:, :)
!! diagonal part of the Hamiltonian, used for preconditioning
COMPLEX(DP) , ALLOCATABLE :: aux2(:, :)
!! temporary storage used in apply_dpot_bands
!
EXTERNAL ch_psi_all, cg_psi
!! functions passed to cgsolve_all
!
! Initialization
!
exclude_hubbard_ = .FALSE.
IF (PRESENT(exclude_hubbard)) exclude_hubbard_ = exclude_hubbard
!
ALLOCATE(h_diag(npwx*npol, nbnd))
ALLOCATE(aux2(npwx*npol, nbnd))
!
all_conv = .TRUE.
tot_num_iter = 0
tot_cg_calls = 0
!
DO ik = 1, nksq
ikk = ikks(ik)
ikq = ikqs(ik)
npw = ngk(ikk)
npwq = ngk(ikq)
!
! Set time-reversed k and k+q points
!
IF (time_reversed) THEN
ikmk = ikmks(ik)
ikmkmq = ikmkmqs(ik)
rsign = -1.0_DP
ELSE
ikmk = ikk
ikmkmq = ikq
rsign = 1.0_DP
ENDIF
!
IF (lsda) current_spin = isk(ikk)
!
! reads unperturbed wavefunctions psi_k in G_space, for all bands
! if q=0, evq is a pointer to evc
!
IF (nksq > 1 .OR. (noncolin .AND. domag)) THEN
IF (lgamma) THEN
CALL get_buffer(evc, lrwfc, iuwfc, ikmk)
ELSE
CALL get_buffer(evc, lrwfc, iuwfc, ikmk)
CALL get_buffer(evq, lrwfc, iuwfc, ikmkmq)
ENDIF
ENDIF
!
! compute beta functions and kinetic energy for k-point ik
! needed by h_psi, called by ch_psi_all, called by cgsolve_all
!
CALL init_us_2(npwq, igk_k(1, ikq), xk(1, ikq), vkb)
CALL g2_kin(ikq)
!
! compute preconditioning matrix h_diag used by cgsolve_all
!
CALL h_prec(ik, evq, h_diag)
!
DO ipert = 1, npert
!
! read P_c^+ x psi_kpoint into dvpsi.
!
nrec = (ipert - 1) * nksq + ik
IF (time_reversed) nrec = nrec + npert * nksq
!
CALL get_buffer(dvpsi, lrdvpsi, iudvpsi, nrec)
!
IF (.NOT. first_iter) THEN
!
! calculates dvscf_q*psi_k in G_space, for all bands, k=kpoint
! dvscf_q from previous iteration (mix_potential)
!
CALL apply_dpot_bands(ik, nbnd_occ(ikk), dvscfins(:, :, ipert), evc, aux2)
dvpsi = dvpsi + aux2
!
! In the case of US pseudopotentials there is an additional
! selfconsist term which comes from the dependence of D on
! V_{eff} on the bare change of the potential
!
IF (time_reversed) THEN
!
! TODO: adddvscf_ph_mag is almost the same as adddvscf, except that it
! uses becp from input, not from USE lrus. Ideally, one should merge the two.
!
CALL adddvscf_ph_mag(ipert, ik, becpt)
ELSE
CALL adddvscf(ipert, ik)
ENDIF
!
! DFPT+U: add to dvpsi the scf part of the response
! Hubbard potential dV_hub
!
IF (lda_plus_u .AND. (.NOT. exclude_hubbard_)) CALL adddvhubscf(ipert, ik)
!
ENDIF
!
! Orthogonalize dvpsi to valence states
!
CALL orthogonalize(dvpsi, evq, ikmk, ikmkmq, dpsi, npwq, .FALSE.)
!
! Initial guess for dpsi
!
IF (first_iter) THEN
!
! At the first iteration dpsi is set to zero
!
dpsi(:, :) = (0.d0,0.d0)
ELSE
!
! starting value for delta_psi is read from iudwf
!
CALL get_buffer(dpsi, lrdwf, iudwf, nrec)
ENDIF
!
! iterative solution of the linear system (H-e)*dpsi=dvpsi
! dvpsi=-P_c+ (dvbare+dvscf)*psi , dvscf fixed.
!
conv_root = .TRUE.
!
! TODO: should nbnd_occ(ikk) be nbnd_occ(ikmk)?
CALL cgsolve_all(ch_psi_all, cg_psi, et(1, ikmk), dvpsi, dpsi, h_diag, &
npwx, npwq, thresh, ik, num_iter, conv_root, anorm, nbnd_occ(ikk), npol)
!
tot_num_iter = tot_num_iter + num_iter
tot_cg_calls = tot_cg_calls + 1
!
IF (.NOT. conv_root) THEN
all_conv = .FALSE.
WRITE( stdout, "(5x, 'kpoint', i4, ' sternheimer_kernel: &
&root not converged, thresh < ', es10.3)") ik, anorm
ENDIF
!
! writes delta_psi on iunit iudwf, k=kpoint,
!
CALL save_buffer(dpsi, lrdwf, iudwf, nrec)
!
! calculates dvscf, sum over k => dvscf_q_ipert
!
IF (noncolin) THEN
CALL incdrhoscf_nc(drhoout(1,1,ipert), wk(ikk), ik, &
dbecsum_nc(1,1,1,1,ipert), dpsi, rsign)
ELSE
CALL incdrhoscf(drhoout(1,current_spin,ipert), wk(ikk), &
ik, dbecsum(1,1,current_spin,ipert), dpsi)
ENDIF
ENDDO ! ipert
ENDDO ! ik
!
CALL mp_sum(tot_num_iter, inter_pool_comm)
CALL mp_sum(tot_cg_calls, inter_pool_comm)
avg_iter = REAL(tot_num_iter, DP) / REAL(tot_cg_calls, DP)
!
!----------------------------------------------------------------------------
END SUBROUTINE sternheimer_kernel
!------------------------------------------------------------------------------
END MODULE response_kernels
!------------------------------------------------------------------------------

13
Makefile
View File

@ -70,7 +70,7 @@ pw : pwlibs
if test -d PW ; then \
( cd PW ; $(MAKE) TLDEPS= all || exit 1) ; fi
cp : bindir libs mods
cp : bindir mods
if test -d CPV ; then \
( cd CPV ; $(MAKE) TLDEPS= all || exit 1) ; fi
@ -112,7 +112,7 @@ gipaw : pwlibs
d3q : phlibs
( cd install ; $(MAKE) -f plugins_makefile $@ || exit 1 )
ld1 : bindir libs mods
ld1 : bindir mods
if test -d atomic ; then \
( cd atomic ; $(MAKE) TLDEPS= all || exit 1 ) ; fi
@ -163,7 +163,7 @@ all : pwall cp ld1 tddfpt hp xspectra gwl
# compile modules, libraries, directory for binaries, etc
###########################################################
pwlibs: bindir libs mods libks_solvers dftd3
pwlibs: bindir mods libks_solvers dftd3
if test -d PW ; then \
( cd PW ; $(MAKE) pw-lib || exit 1) ; fi
@ -186,7 +186,7 @@ pw4gwwlib : phlibs
mods : libfox libutil libla libfft libupf libmbd librxc
( cd Modules ; $(MAKE) TLDEPS= all || exit 1 )
libks_solvers : libs libutil libla
libks_solvers : libutil libla
( cd KS_Solvers ; $(MAKE) TLDEPS= all || exit 1 )
libla : liblapack libutil libcuda
@ -204,9 +204,6 @@ libutil :
libupf : libutil libcuda
( cd upflib ; $(MAKE) TLDEPS= all || exit 1 )
libs :
( cd clib ; $(MAKE) TLDEPS= all || exit 1 )
lrmods : mods pwlibs
( cd LR_Modules ; $(MAKE) TLDEPS= all || exit 1 )
@ -287,7 +284,7 @@ clean :
for dir in \
CPV LAXlib FFTXlib XClib UtilXlib upflib Modules PP PW EPW KS_Solvers \
NEB ACFDT COUPLE GWW XSpectra PWCOND dft-d3 \
atomic clib LR_Modules upflib \
atomic LR_Modules upflib \
dev-tools extlibs Environ TDDFPT PHonon HP GWW Doc GUI \
QEHeat \
; do \

16
Modules/CMakeLists.txt
View File

@ -93,7 +93,8 @@ set(src_modules
wyckoff.f90
wypos.f90
zvscal.f90
# list of subroutines and functions (not modules) previously found in flib/
wave_gauge.f90
# subroutines and functions (not modules) previously found in flib
atom_weight.f90
capital.f90
cryst_to_car.f90
@ -129,10 +130,19 @@ set(src_modules
random_numbers_gpu.f90
cuda_subroutines.f90)
qe_enable_cuda_fortran("${src_modules}")
qe_add_library(qe_modules ${src_modules})
# subroutines and functions (not modules) previously found in clib
set(src_modules_c
customize_signals.c
qmmm_aux.c
sockets.c
stack.c)
qe_add_library(qe_modules_c ${src_modules_c})
target_link_libraries(qe_modules
PRIVATE
qe_modules_c
qe_openmp_fortran
qe_fftx
qe_lax
@ -153,4 +163,4 @@ if(QE_ENABLE_CUDA)
endif()
###########################################################
qe_install_targets(qe_modules)
qe_install_targets(qe_modules qe_modules_c)

11
Modules/Makefile
View File

@ -102,9 +102,10 @@ fox_init_module.o \
xsf.o \
wyckoff.o \
wypos.o \
zvscal.o
zvscal.o \
wave_gauge.o
# list of subroutines and functions (not modules) previously found in flib/
# list of subroutines and functions (not modules) previously found in flib/clib
OBJS = \
atom_weight.o \
@ -133,7 +134,11 @@ volume.o \
wgauss.o \
w0gauss.o \
w1gauss.o \
deviatoric.o
deviatoric.o \
customize_signals.o \
qmmm_aux.o \
sockets.o \
stack.o
# GPU versions of modules
MODULES += \

0
clib/customize_signals.c → Modules/customize_signals.c
View File

3
Modules/make.depend
View File

@ -395,6 +395,9 @@ wannier_gw.o : io_global.o
wannier_gw.o : kind.o
wannier_gw.o : recvec.o
wannier_new.o : kind.o
wave_gauge.o : ../UtilXlib/mp.o
wave_gauge.o : kind.o
wave_gauge.o : mp_pools.o
wavefunctions.o : kind.o
wavefunctions_gpu.o : wavefunctions.o
wgauss.o : constants.o

3
Modules/qexsd.f90
View File

@ -585,6 +585,7 @@ SUBROUTINE qexsd_init_clocks (timing_, total_clock, partial_clocks)
match = clock_label(1:nclock) == TRIM(partial_clocks(ipar))
IF ( ANY (match)) THEN
nc = get_index(.TRUE., match)
IF (nc == ic) CYCLE
t = get_cpu_and_wall(nc)
CALL qes_init(partial_(ipar), "partial", TRIM(clock_label(nc)), t(1), t(2),&
called(nc))
@ -592,7 +593,7 @@ SUBROUTINE qexsd_init_clocks (timing_, total_clock, partial_clocks)
CALL qes_init (partial_(ipar), "partial", "not_found", -1.d0, -1.d0, 0)
partial_(ipar)%lwrite=.FALSE.
END IF
END DO
END DO
END IF
CALL qes_init( timing_, "timing_info", total_, partial_)
CALL qes_reset ( total_)

0
clib/qmmm_aux.c → Modules/qmmm_aux.c
View File

0
clib/sockets.c → Modules/sockets.c
View File

0
clib/stack.c → Modules/stack.c
View File

155
Modules/wave_gauge.f90 Normal file
View File

@ -0,0 +1,155 @@
!
! Copyright (C) 2001-2021 Quantum ESPRESSO group
! This file is distributed under the terms of the
! GNU General Public License. See the file `License'
! in the root directory of the present distribution,
! or http://www.gnu.org/copyleft/gpl.txt .
!
!
!-----------------------------------------------------------------------
MODULE wave_gauge
! this module contains routines that are used to compute the numerical
! derivative of the wavefunctions fixing a particular gauge
USE kinds, ONLY: DP
implicit none
contains
subroutine compute_dot_evc_parallel_gauge(t_minus, t_zero, t_plus, dot_evc, nbnd, npw, npwx, gstart)
! compute the numerical derivative using the wavefunctions computed at t-dt/2, t, t+dt/2
! in the parallel transport gauge
integer, intent(in) :: nbnd, npw, npwx, gstart
COMPLEX(DP), intent(inout) :: dot_evc(:, :)
complex(dp), intent(in) :: t_minus(:,:), t_zero(:,:), t_plus(:,:)
dot_evc=(0.0_dp, 0.0_dp)
call project_parallel_gauge(t_minus, t_zero, t_plus, dot_evc, nbnd, npw, npwx, gstart, 1.0_dp,&
.true., .true.)
end subroutine
subroutine project_parallel_gauge_2(t_minus, t_zero, t_zero_proj , nbnd, npw, npwx, gstart)
! project t_zero over the manifold of t_minus in the parallel transport gauge.
! Do not project over the not occupied manifold
integer, intent(in) :: nbnd, npw, npwx, gstart
COMPLEX(DP), intent(inout) :: t_zero_proj(:,:)
complex(dp), intent(in) :: t_minus(:, :), t_zero(:,:)
complex(dp), allocatable :: dummy(:,:)
t_zero_proj=0.0_dp
call project_parallel_gauge(t_zero, t_minus, dummy, t_zero_proj, nbnd, npw, npwx,&
gstart, -1.0_dp, .false., .false. )
end subroutine
subroutine project_parallel_gauge(t_minus, t_zero, t_plus, dot_evc, nbnd, npw, npwx, gstart, &
factor, use_t_plus, project_conduction)
!! implemented only in the plain dft norm conserving pseudopotential case
!! let P be the projector over the occupied manifold space:
!! P = \sum_v |c_v><c_v|
!! then (no surprise)
!! |c> = P |c>
!! doing a time derivative:
!! \dot |c> = \dot P |c> + P \dot |c>
!! the parallel transport gauge fixes P \dot |c> to zero (the derivative
!! has no components over the occupied space). \dot P can be computed numerically
!! using wavefunctions that possibly have different gauges.
!! Then, to make sure that there are no components over the occupied space,
!! a projector over the conduction band is added if project_conduction is .true.:
!!
!! (1-P) \dot P |c>
!!
!! t_minus, t_zero, t_plus are the wavefunctions. t_minus and t_plus are separated by a time dt
!! t_zero is computed in the middle. The gauge of the result is the same of t_zero.
!! This works also if you put t_zero = t_minus or t_zero = t_plus (be careful to the gauge!).
!! This routine adds \dot |c>*dt*factor (the derivative multiplied by dt multiplied by the specified factor)
!! to the output array dot_evc with the gauge of t_zero. Gauges of t_minus and t_plus should not enter the
!! result and can be arbitrary. If use_t_plus is false, t_plus = t_zero is assumed and some cpu time
!! is saved. If project_conduction is true the (1-P) projector is computed, and this involves an additional
!! matrix-matrix product with dimensions nbnd x nbnd
use mp, only: mp_sum
USE mp_pools, ONLY: intra_pool_comm
integer, intent(in) :: nbnd, npw, npwx, gstart
COMPLEX(DP), intent(inout) :: dot_evc(:, :)
complex(dp), intent(in) :: t_minus(:,:), t_zero(:,:), t_plus(:,:)
real(dp), intent(in) :: factor
logical, intent(in) :: use_t_plus, project_conduction
real(DP), allocatable :: sa(:,:), ssa(:,:), sb(:,:), ssb(:,:)
complex(dp) :: tmp
integer :: ibnd, jbnd, ig
allocate(sb(nbnd, nbnd))
! \dot P |c> =
! computed sa_ij = < t_plus_i, t_zero_j >, sb_ij = < t_minus_i, t_zero_j >
! remove contribution at G=0 due to c(-g)*=c(g), and only half vector is stored
! (gamma trick)
! sa is the identity if 2 points are used (t_plus==t_zero)
if (use_t_plus) then
allocate(sa(nbnd, nbnd))
call dgemm('T', 'N', nbnd, nbnd, 2*npw, 2.d0, t_plus, 2*npwx, t_zero, 2*npwx, 0.d0, sa, nbnd)
end if
call dgemm('T', 'N', nbnd, nbnd, 2*npw, 2.d0, t_minus, 2*npwx, t_zero, 2*npwx, 0.d0, sb, nbnd)
if (gstart == 2) then
do ibnd = 1, nbnd
do jbnd = 1, nbnd
if (use_t_plus) &
sa(ibnd, jbnd) = sa(ibnd, jbnd) - dble(conjg(t_plus(1, ibnd))*t_zero(1, jbnd))
sb(ibnd, jbnd) = sb(ibnd, jbnd) - dble(conjg(t_minus(1, ibnd))*t_zero(1, jbnd))
end do
end do
end if
call mp_sum(sb, intra_pool_comm)
if (use_t_plus) &
call mp_sum(sa, intra_pool_comm)
! compute scalar products that appear due to the ( 1 - P ) projector over
! the empty bands, that are
! ssa = sa ^ 2 and ssb = sb ^ 2. The gauge cancels out here
if (project_conduction) then
allocate(ssb(nbnd, nbnd))
if (use_t_plus) then
allocate(ssa(nbnd, nbnd))
call dgemm('T', 'N', nbnd, nbnd, nbnd, 1.d0, sa, nbnd, sa, nbnd, 0.d0, ssa, nbnd)
end if
call dgemm('T', 'N', nbnd, nbnd, nbnd, 1.d0, sb, nbnd, sb, nbnd, 0.d0, ssb, nbnd)
end if
! compute final projection
do ibnd = 1, nbnd
do jbnd = 1, nbnd
do ig = 1, npw
tmp=0.0_dp
tmp = tmp - t_minus(ig, jbnd)*sb(jbnd, ibnd)
if (use_t_plus) then
tmp = tmp + t_plus(ig, jbnd)*sa(jbnd, ibnd)
endif
if (project_conduction) then
if (use_t_plus)&
tmp = tmp - t_zero(ig, jbnd)*ssa(ibnd, jbnd)
tmp = tmp + t_zero(ig, jbnd)*ssb(ibnd, jbnd)
endif
dot_evc(ig, ibnd) = dot_evc(ig, ibnd) + factor*tmp
end do
end do
end do
if (project_conduction) then
deallocate(ssb)
if (use_t_plus)&
deallocate(ssa)
end if
if (use_t_plus)&
deallocate(sa)
deallocate(sb)
end subroutine
end module

14
NEB/src/Makefile
View File

@ -38,22 +38,20 @@ path_variables.o \
set_defaults.o \
stop_run_path.o
PWOBJS= ../../PW/src/libpw.a
QEMODS = ../../Modules/libqemod.a ../../KS_Solvers/libks_solvers.a ../../upflib/libupf.a \
../../FFTXlib/libqefft.a ../../LAXlib/libqela.a ../../UtilXlib/libutil.a ../../dft-d3/libdftd3qe.a \
../../XClib/xc_lib.a
PWOBJS= ../../PW/src/libpw.a ../../KS_Solvers/libks_solvers.a ../../dft-d3/libdftd3qe.a
QEMODS = $(BASEMODS)
TLDEPS=pwlibs
all : tldeps neb.x path_interpolation.x
neb.x : $(NEBOBJS) libneb.a $(LIBOBJS) $(PWOBJS) $(QEMODS)
neb.x : $(NEBOBJS) libneb.a $(PWOBJS) $(QEMODS)
$(LD) $(LDFLAGS) -o $@ \
$(NEBOBJS) libneb.a $(PWOBJS) $(QEMODS) $(LIBOBJS) $(QELIBS)
$(NEBOBJS) libneb.a $(PWOBJS) $(QEMODS) $(QELIBS)
- ( cd ../../bin; ln -fs ../NEB/src/$@ . )
path_interpolation.x : path_interpolation.o libneb.a $(PWOBJS) $(QEMODS) $(LIBOBJS)
path_interpolation.x : path_interpolation.o libneb.a $(PWOBJS) $(QEMODS)
$(LD) $(LDFLAGS) -o $@ \
path_interpolation.o libneb.a $(PWOBJS) $(QEMODS) $(LIBOBJS) $(QELIBS)
path_interpolation.o libneb.a $(PWOBJS) $(QEMODS) $(QELIBS)
- ( cd ../../bin ; ln -fs ../NEB/src/$@ . )
libneb.a : $(NEBLIBS)

2
PHonon/CMakeLists.txt
View File

@ -10,7 +10,6 @@ set(src_ph
PH/addusddens.f90
PH/addusddense.f90
PH/addusdynmat.f90
PH/adddvhubscf.f90
PH/allocate_part.f90
PH/allocate_pert.f90
PH/allocate_phq.f90
@ -178,7 +177,6 @@ set(src_ph
PH/dhdrhopsi.f90
PH/set_dvscf.f90
PH/non_scf_ph.f90
PH/adddvscf_ph_mag.f90
PH/apply_trev.f90
PH/c_bands_ph.f90
PH/dvscf_interpolate.f90

8
PHonon/Doc/INPUT_PH.def
View File

@ -166,9 +166,11 @@ input_description -distribution {Quantum ESPRESSO} -package PHonon -program ph.x
var trans -type LOGICAL {
default { .true. }
info {
If .true. the phonons are computed.
If @ref trans .and. @ref epsil are .true. effective charges are
calculated.
If .false. the phonons are not computed.
If @ref trans .and. @ref epsil are both .true.,
the effective charges are calculated.
If @ref ldisp is .true., trans=.false. is overridden
(except for the case of electron-phonon calculations)
}
}

13
PHonon/FD/Makefile
View File

@ -11,11 +11,8 @@ FDOBJS = \
stop_pp.o
QEMODS = ../../Modules/libqemod.a ../../FFTXlib/libqefft.a ../../upflib/libupf.a \
../../KS_Solvers/libks_solvers.a ../../XClib/xc_lib.a \
../../LAXlib/libqela.a ../../UtilXlib/libutil.a ../../dft-d3/libdftd3qe.a
PWOBJS = ../../PW/src/libpw.a
PWOBJS = ../../PW/src/libpw.a ../../KS_Solvers/libks_solvers.a ../../dft-d3/libdftd3qe.a
QEMODS = $(BASEMODS)
TLDEPS= pwlibs
@ -23,17 +20,17 @@ all : tldeps fd.x fd_ifc.x fd_ef.x
fd_ifc.x : fd_ifc.o libfd.a $(FDOBJS) $(PWOBJS) $(QEMODS)
$(LD) $(LDFLAGS) -o $@ \
fd_ifc.o libfd.a $(PWOBJS) $(QEMODS) $(LIBOBJS) $(QELIBS)
fd_ifc.o libfd.a $(PWOBJS) $(QEMODS) $(QELIBS)
- ( cd ../../bin ; ln -fs ../PHonon/FD/$@ . )
fd.x : fd.o libfd.a $(FDOBJS) $(PWOBJS) $(QEMODS)
$(LD) $(LDFLAGS) -o $@ \
fd.o libfd.a $(PWOBJS) $(QEMODS) $(LIBOBJS) $(QELIBS)
fd.o libfd.a $(PWOBJS) $(QEMODS) $(QELIBS)
- ( cd ../../bin ; ln -fs ../PHonon/FD/$@ . )
fd_ef.x : fd_ef.o libfd.a $(FDOBJS) $(PWOBJS) $(QEMODS)
$(LD) $(LDFLAGS) -o $@ \
fd_ef.o libfd.a $(PWOBJS) $(QEMODS) $(LIBOBJS) $(QELIBS)
fd_ef.o libfd.a $(PWOBJS) $(QEMODS) $(QELIBS)
- ( cd ../../bin ; ln -fs ../PHonon/FD/$@ . )
tldeps :
if test -n "$(TLDEPS)" ; then \

9
PHonon/Gamma/Makefile
View File

@ -34,17 +34,16 @@ solve_e.o \
solve_ph.o \
writedyn.o
PWOBJS = ../../PW/src/libpw.a
QEMODS = ../../Modules/libqemod.a ../../KS_Solvers/libks_solvers.a ../../upflib/libupf.a ../../XClib/xc_lib.a \
../../FFTXlib/libqefft.a ../../LAXlib/libqela.a ../../UtilXlib/libutil.a ../../dft-d3/libdftd3qe.a ../../LR_Modules/liblrmod.a
PWOBJS = ../../PW/src/libpw.a ../../KS_Solvers/libks_solvers.a ../../dft-d3/libdftd3qe.a
QEMODS = ../../LR_Modules/liblrmod.a $(BASEMODS)
TLDEPS= pwlibs
all : tldeps phcg.x
phcg.x : phcg.o libphcg.a $(PWOBJS) $(QEMODS) $(LIBOBJS)
phcg.x : phcg.o libphcg.a $(PWOBJS) $(QEMODS)
$(LD) $(LDFLAGS) -o $@ \
phcg.o libphcg.a $(PWOBJS) $(QEMODS) $(LIBOBJS) $(QELIBS)
phcg.o libphcg.a $(PWOBJS) $(QEMODS) $(QELIBS)
- ( cd ../../bin ; ln -fs ../PHonon/Gamma/phcg.x . )
tldeps :
if test -n "$(TLDEPS)" ; then \

51
PHonon/PH/Makefile
View File

@ -14,13 +14,11 @@ add_for_charges.o \
add_zstar_ue.o \
add_zstar_ue_us.o \
addcore.o \
adddvscf_ph_mag.o \
addnlcc.o \
addnlcc_zstar_eu_us.o \
addusddens.o \
addusddense.o \
addusdynmat.o \
adddvhubscf.o \
allocate_part.o \
allocate_pert.o \
allocate_phq.o \
@ -198,9 +196,8 @@ PHAUXOBJS = \
write_eigenvectors.o
LRMODS = ../../LR_Modules/liblrmod.a
PWOBJS = ../../PW/src/libpw.a
QEMODS = ../../Modules/libqemod.a ../../KS_Solvers/libks_solvers.a ../../upflib/libupf.a ../../XClib/xc_lib.a \
../../FFTXlib/libqefft.a ../../LAXlib/libqela.a ../../UtilXlib/libutil.a ../../dft-d3/libdftd3qe.a
PWOBJS = ../../PW/src/libpw.a ../../KS_Solvers/libks_solvers.a ../../dft-d3/libdftd3qe.a
QEMODS = $(BASEMODS)
TLDEPS= phlibs
all : tldeps libs-ph ph.x dynmat.x matdyn.x q2r.x lambda.x fqha.x q2qstar.x \
@ -208,63 +205,63 @@ all : tldeps libs-ph ph.x dynmat.x matdyn.x q2r.x lambda.x fqha.x q2qstar.x \
libs-ph : libph.a libphaux.a
ph.x : phonon.o libph.a $(PWOBJS) $(LRMODS) $(QEMODS) $(LIBOBJS)
ph.x : phonon.o libph.a $(PWOBJS) $(LRMODS) $(QEMODS)
$(LD) $(LDFLAGS) -o $@ phonon.o libph.a \
$(LRMODS) $(PWOBJS) $(QEMODS) $(LIBOBJS) $(QELIBS)
$(LRMODS) $(PWOBJS) $(QEMODS) $(QELIBS)
- ( cd ../../bin ; ln -fs ../PHonon/PH/ph.x . )
dynmat.x : dynmat.o libphaux.a libph.a $(PWOBJS) $(LRMODS) $(QEMODS) $(LIBOBJS)
dynmat.x : dynmat.o libphaux.a libph.a $(PWOBJS) $(LRMODS) $(QEMODS)
$(LD) $(LDFLAGS) -o $@ dynmat.o libphaux.a libph.a \
$(LRMODS) $(PWOBJS) $(QEMODS) $(LIBOBJS) $(QELIBS)
$(LRMODS) $(PWOBJS) $(QEMODS) $(QELIBS)
- ( cd ../../bin ; ln -fs ../PHonon/PH/$@ . )
matdyn.x : matdyn.o libphaux.a libph.a $(PWOBJS) $(LRMODS) $(QEMODS) $(LIBOBJS)
matdyn.x : matdyn.o libphaux.a libph.a $(PWOBJS) $(LRMODS) $(QEMODS)
$(LD) $(LDFLAGS) -o $@ matdyn.o libphaux.a libph.a \
$(LRMODS) $(PWOBJS) $(QEMODS) $(LIBOBJS) $(QELIBS)
$(LRMODS) $(PWOBJS) $(QEMODS) $(QELIBS)
- ( cd ../../bin ; ln -fs ../PHonon/PH/$@ . )
q2r.x : q2r.o libphaux.a libph.a $(PWOBJS) $(LRMODS) $(QEMODS) $(LIBOBJS)
q2r.x : q2r.o libphaux.a libph.a $(PWOBJS) $(LRMODS) $(QEMODS)
$(LD) $(LDFLAGS) -o $@ q2r.o libphaux.a libph.a \
$(LRMODS) $(PWOBJS) $(QEMODS) $(LIBOBJS) $(QELIBS)
$(LRMODS) $(PWOBJS) $(QEMODS) $(QELIBS)
- ( cd ../../bin ; ln -fs ../PHonon/PH/$@ . )
q2qstar.x : q2qstar.o libph.a $(PWOBJS) $(LRMODS) $(QEMODS) $(LIBOBJS)
q2qstar.x : q2qstar.o libph.a $(PWOBJS) $(LRMODS) $(QEMODS)
$(LD) $(LDFLAGS) -o $@ q2qstar.o libph.a \
$(LRMODS) $(PWOBJS) $(QEMODS) $(LIBOBJS) $(QELIBS)
$(LRMODS) $(PWOBJS) $(QEMODS) $(QELIBS)
- ( cd ../../bin ; ln -fs ../PHonon/PH/$@ . )
lambda.x : lambda.o $(PWOBJS) $(LRMODS) $(QEMODS) $(LIBOBJS)
lambda.x : lambda.o $(PWOBJS) $(LRMODS) $(QEMODS)
$(LD) $(LDFLAGS) -o $@ lambda.o \
$(LRMODS) $(PWOBJS) $(QEMODS) $(LIBOBJS) $(QELIBS)
$(LRMODS) $(PWOBJS) $(QEMODS) $(QELIBS)
- ( cd ../../bin ; ln -fs ../PHonon/PH/$@ . )
alpha2f.x : alpha2f.o libph.a $(PWOBJS) $(LRMODS) $(QEMODS) $(LIBOBJS)
alpha2f.x : alpha2f.o libph.a $(PWOBJS) $(LRMODS) $(QEMODS)
$(LD) $(LDFLAGS) -o $@ alpha2f.o libph.a \
$(LRMODS) $(PWOBJS) $(QEMODS) $(LIBOBJS) $(QELIBS)
$(LRMODS) $(PWOBJS) $(QEMODS) $(QELIBS)
- ( cd ../../bin ; ln -fs ../PHonon/PH/$@ . )
epa.x : epa.o libph.a $(PWOBJS) $(LRMODS) $(QEMODS) $(LIBOBJS)
epa.x : epa.o libph.a $(PWOBJS) $(LRMODS) $(QEMODS)
$(LD) $(LDFLAGS) -o $@ epa.o libph.a \
$(LRMODS) $(PWOBJS) $(QEMODS) $(LIBOBJS) $(QELIBS)
$(LRMODS) $(PWOBJS) $(QEMODS) $(QELIBS)
- ( cd ../../bin ; ln -fs ../PHonon/PH/$@ . )
dvscf_q2r.x : dvscf_q2r.o libphaux.a libph.a $(PWOBJS) $(LRMODS) $(QEMODS) $(LIBOBJS)
dvscf_q2r.x : dvscf_q2r.o libphaux.a libph.a $(PWOBJS) $(LRMODS) $(QEMODS)
$(LD) $(LDFLAGS) -o $@ dvscf_q2r.o libphaux.a libph.a \
$(LRMODS) $(PWOBJS) $(QEMODS) $(LIBOBJS) $(QELIBS)
$(LRMODS) $(PWOBJS) $(QEMODS) $(QELIBS)
- ( cd ../../bin ; ln -fs ../PHonon/PH/$@ . )
#fqha.o :
# $(MPIF90) $(FFLAGS_NOOPT) -c fqha.f90
fqha.x : fqha.o $(PWOBJS) $(QEMODS) $(LIBOBJS)
fqha.x : fqha.o $(PWOBJS) $(QEMODS)
$(LD) $(LDFLAGS) -o $@ fqha.o \
$(PWOBJS) $(QEMODS) $(LIBOBJS) $(QELIBS)
$(PWOBJS) $(QEMODS) $(QELIBS)
- ( cd ../../bin ; ln -fs ../PHonon/PH/$@ . )
postahc.x : postahc.o $(PWOBJS) $(QEMODS) $(LIBOBJS)
postahc.x : postahc.o $(PWOBJS) $(QEMODS)
$(LD) $(LDFLAGS) -o $@ postahc.o \
$(PWOBJS) $(QEMODS) $(LIBOBJS) $(QELIBS)
$(PWOBJS) $(QEMODS) $(QELIBS)
- ( cd ../../bin ; ln -fs ../PHonon/PH/$@ . )
tldeps :

3
PHonon/PH/add_zstar_ue.f90
View File

@ -23,8 +23,7 @@ subroutine add_zstar_ue (imode0, npe)
USE qpoint, ONLY: nksq, ikks
USE eqv, ONLY: dpsi, dvpsi
USE efield_mod, ONLY: zstarue0_rec
USE units_ph, ONLY : iudwf, lrdwf
USE units_lr, ONLY : iuwfc, lrwfc
USE units_lr, ONLY : iuwfc, lrwfc, iudwf, lrdwf
USE control_lr, ONLY : nbnd_occ
USE uspp_init, ONLY : init_us_2

4
PHonon/PH/allocate_phq.f90
View File

@ -46,8 +46,8 @@ subroutine allocate_phq
USE ldaU, ONLY : lda_plus_u, Hubbard_lmax, nwfcU
USE ldaU_ph, ONLY : dnsbare, dnsorth, dnsbare_all_modes, wfcatomk, &
dwfcatomk, sdwfcatomk, wfcatomkpq, dwfcatomkpq, &
swfcatomk, swfcatomkpq, sdwfcatomkpq, dvkb, vkbkpq, &
dvkbkpq
sdwfcatomkpq, dvkb, vkbkpq, dvkbkpq
USE ldaU_lr, ONLY : swfcatomk, swfcatomkpq
USE qpoint_aux, ONLY : becpt, alphapt
IMPLICIT NONE

3
PHonon/PH/chi_test.f90
View File

@ -18,8 +18,7 @@ subroutine chi_test (dvscfs, chif, ik, depsi, auxr, auxg)
USE wvfct, ONLY : npwx, nbnd
USE fft_base, ONLY : dffts
use ramanm, ONLY : lrd2w, iud2w, jab
USE units_ph, ONLY : lrdwf, iudwf
USE units_lr, ONLY : iuwfc
USE units_lr, ONLY : lrdwf, iudwf
USE buffers, ONLY : get_buffer
USE qpoint, ONLY : npwq, nksq
USE eqv, ONLY : dpsi, dvpsi

4
PHonon/PH/close_phq.f90
View File

@ -19,10 +19,10 @@ SUBROUTINE close_phq( flag )
USE io_files, ONLY : iunhub, iunhub_noS
USE buffers, ONLY : close_buffer
USE uspp, ONLY : okvan
USE units_ph, ONLY : iudwf, iubar, iudrhous, iuebar, iudrho, &
USE units_ph, ONLY : iubar, iudrhous, iuebar, iudrho, &
iudvscf, iucom, iudvkb3, iuint3paw, iudyn, &
iundnsscf, iudvpsi, iugauge
USE units_lr, ONLY : iuwfc, iuatwfc, iuatswfc
USE units_lr, ONLY : iuwfc, iuatwfc, iuatswfc, iudwf
USE control_ph, ONLY : zue, epsil, only_wfc
USE recover_mod, ONLY : clean_recover
USE output, ONLY : fildrho, fildvscf

6
PHonon/PH/d2nsq_bare.f90
View File

@ -1078,8 +1078,9 @@ SUBROUTINE d2nsq_bare_k (ik, icart, jcart, na, nap, nah, &
USE units_lr, ONLY : iuwfc, lrwfc
USE ions_base, ONLY : nat, ityp, ntyp => nsp
USE klist, ONLY : xk, ngk, igk_k
USE ldaU_ph, ONLY : wfcatomk, swfcatomk, wfcatomkpq, dwfcatomk, dwfcatomkpq, &
USE ldaU_ph, ONLY : wfcatomk, wfcatomkpq, dwfcatomk, dwfcatomkpq, &
dvkb, vkbkpq, dvkbkpq, proj1, d2ns_type
USE ldaU_lr, ONLY : swfcatomk
USE wvfct, ONLY : npwx, nbnd, wg
USE uspp, ONLY : vkb, nkb
USE qpoint, ONLY : nksq, ikks, ikqs
@ -1258,8 +1259,9 @@ SUBROUTINE d2nsq_bare_k_diag (ik, icart, jcart, na, nap, nah, &
USE units_lr, ONLY : iuwfc, lrwfc
USE ions_base, ONLY : nat, ityp, ntyp => nsp
USE klist, ONLY : xk, ngk, igk_k
USE ldaU_ph, ONLY : wfcatomk, swfcatomk, wfcatomkpq, dwfcatomk, dwfcatomkpq, &
USE ldaU_ph, ONLY : wfcatomk, wfcatomkpq, dwfcatomk, dwfcatomkpq, &
dvkb, vkbkpq, dvkbkpq, proj1, d2ns_type
USE ldaU_lr, ONLY : swfcatomk
USE wvfct, ONLY : npwx, nbnd, wg
USE uspp, ONLY : vkb, nkb
USE qpoint, ONLY : nksq, ikks, ikqs

5
PHonon/PH/deallocate_phq.f90
View File

@ -39,9 +39,10 @@ subroutine deallocate_phq
USE control_lr, ONLY : lgamma, nbnd_occ
USE ldaU, ONLY : lda_plus_u
USE ldaU_ph, ONLY : dnsbare_all_modes, dnsorth_cart, dnsorth, dnsbare, &
wfcatomk, swfcatomk, dwfcatomk, sdwfcatomk, &
wfcatomkpq, dwfcatomkpq, swfcatomkpq, sdwfcatomkpq, &
wfcatomk, dwfcatomk, sdwfcatomk, &
wfcatomkpq, dwfcatomkpq, sdwfcatomkpq, &
dvkb, vkbkpq, dvkbkpq
USE ldaU_lr, ONLY : swfcatomk, swfcatomkpq
USE qpoint_aux, ONLY : ikmks, ikmkmqs, becpt, alphapt
USE becmod, ONLY : deallocate_bec_type
USE nc_mag_aux, ONLY : int1_nc_save, deeq_nc_save

3
PHonon/PH/dhdrhopsi.f90
View File

@ -47,8 +47,7 @@ subroutine dhdrhopsi
USE becmod, ONLY : calbec, bec_type, allocate_bec_type, &
deallocate_bec_type, beccopy
use ramanm, ONLY : lrchf, iuchf, lrd2w, iud2w, jab, dek, eth_ns
USE units_ph, ONLY : lrdwf, iudwf
USE units_lr, ONLY : iuwfc, lrwfc
USE units_lr, ONLY : iuwfc, lrwfc, lrdwf, iudwf
USE lrus, ONLY : becp1
USE eqv, ONLY : dpsi, dvpsi
USE qpoint, ONLY : nksq

3
PHonon/PH/dielec.f90
View File

@ -22,7 +22,8 @@ subroutine dielec()
USE buffers, ONLY : get_buffer
USE noncollin_module, ONLY : npol
USE efield_mod, ONLY : epsilon
USE units_ph, ONLY : lrdwf, iudwf, lrebar, iuebar
USE units_ph, ONLY : lrebar, iuebar
USE units_lr, ONLY : lrdwf, iudwf
USE eqv, ONLY : dpsi, dvpsi
USE qpoint, ONLY : nksq, ikks
USE ph_restart, ONLY : ph_writefile

3
PHonon/PH/dnsq_bare.f90
View File

@ -29,9 +29,10 @@ SUBROUTINE dnsq_bare
USE ions_base, ONLY : nat, ityp, ntyp => nsp
USE klist, ONLY : xk, ngk, igk_k
USE ldaU, ONLY : Hubbard_lmax, Hubbard_l, offsetU, is_hubbard, nwfcU
USE ldaU_ph, ONLY : wfcatomk, wfcatomkpq, swfcatomk, swfcatomkpq, dwfcatomkpq, &
USE ldaU_ph, ONLY : wfcatomk, wfcatomkpq, dwfcatomkpq, &
sdwfcatomk, sdwfcatomkpq, dvkb, vkbkpq, dvkbkpq, &
dnsbare, dnsbare_all_modes, proj1, proj2, read_dns_bare
USE ldaU_lr, ONLY : swfcatomk, swfcatomkpq
USE wvfct, ONLY : npwx, wg, nbnd
USE uspp, ONLY : vkb, nkb, ofsbeta
USE qpoint, ONLY : nksq, ikks, ikqs

3
PHonon/PH/dnsq_orth.f90
View File

@ -30,9 +30,10 @@ SUBROUTINE dnsq_orth()
USE units_lr, ONLY : iuwfc, lrwfc
USE ions_base, ONLY : nat, ityp, ntyp => nsp
USE ldaU, ONLY : Hubbard_lmax, Hubbard_l, is_hubbard, offsetU, nwfcU
USE ldaU_ph, ONLY : swfcatomk, swfcatomkpq, dvkb, vkbkpq, dvkbkpq, &
USE ldaU_ph, ONLY : dvkb, vkbkpq, dvkbkpq, &
proj1, proj2, dnsorth_cart, &
read_dns_bare, dnsorth
USE ldaU_lr, ONLY : swfcatomk, swfcatomkpq
USE klist, ONLY : xk, wk, ngk, igk_k
USE wvfct, ONLY : npwx, wg, nbnd
USE qpoint, ONLY : nksq, ikks, ikqs

6
PHonon/PH/dnsq_scf.f90
View File

@ -25,11 +25,11 @@ SUBROUTINE dnsq_scf (npe, lmetq0, imode0, irr, lflag)
!
USE kinds, ONLY : DP
USE io_files, ONLY : nwordwfcU
USE units_ph, ONLY : lrdwf, iudwf
USE units_lr, ONLY : iuwfc, lrwfc
USE units_lr, ONLY : iuwfc, lrwfc, iudwf, lrdwf
USE ions_base, ONLY : nat, ityp, ntyp => nsp
USE ldaU, ONLY : Hubbard_lmax, Hubbard_l, is_hubbard, offsetU, nwfcU
USE ldaU_ph, ONLY : swfcatomk, swfcatomkpq, proj1, proj2, dnsscf, dnsorth
USE ldaU_ph, ONLY : proj1, proj2, dnsorth
USE ldaU_lr, ONLY : swfcatomk, swfcatomkpq, dnsscf
USE klist, ONLY : xk, wk, degauss, ngauss, ngk
USE wvfct, ONLY : npwx, wg, nbnd, et
USE qpoint, ONLY : nksq, ikks, ikqs

2
PHonon/PH/drhodv.f90
View File

@ -39,7 +39,7 @@ subroutine drhodv (nu_i0, nper, drhoscf)
USE dynmat, ONLY : dyn, dyn_rec
USE modes, ONLY : u
USE units_ph, ONLY : lrdwf, iudwf
USE units_lr, ONLY : lrdwf, iudwf
USE eqv, ONLY : dpsi
USE qpoint, ONLY : nksq, ikks, ikqs

4
PHonon/PH/dvpsi_e2.f90
View File

@ -25,8 +25,8 @@ subroutine dvpsi_e2
USE fft_interfaces, ONLY : fft_interpolate
USE scf, ONLY : rho
USE qpoint, ONLY : nksq
USE units_ph, ONLY : lrdrho, iudrho, lrdwf, iudwf
USE units_lr, ONLY : iuwfc, lrwfc
USE units_ph, ONLY : lrdrho, iudrho
USE units_lr, ONLY : iuwfc, lrwfc, lrdwf, iudwf
USE control_lr, ONLY : nbnd_occ
USE ramanm, ONLY : lrba2, iuba2, lrchf, iuchf, a1j, a2j
USE mp_pools, ONLY : my_pool_id, inter_pool_comm

5
PHonon/PH/dvqhub_barepsi_us.f90
View File

@ -43,9 +43,10 @@ SUBROUTINE dvqhub_barepsi_us (ik, uact)
USE ions_base, ONLY : nat, ityp, ntyp => nsp
USE klist, ONLY : xk, ngk, igk_k
USE ldaU, ONLY : U_projection, Hubbard_l, is_hubbard, Hubbard_J0, offsetU, nwfcU
USE ldaU_ph, ONLY : wfcatomk, wfcatomkpq, swfcatomk, swfcatomkpq, dwfcatomkpq, &
USE ldaU_ph, ONLY : wfcatomk, wfcatomkpq, dwfcatomkpq, &
sdwfcatomk, sdwfcatomkpq, dvkb, vkbkpq, dvkbkpq, &
proj1, proj2, dnsbare, effU
proj1, proj2, dnsbare
USE ldaU_lr, ONLY : effU, swfcatomk, swfcatomkpq
USE wvfct, ONLY : npwx, nbnd
USE uspp, ONLY : vkb, nkb, ofsbeta
USE qpoint, ONLY : nksq, ikks, ikqs

5
PHonon/PH/dvqhub_barepsi_us2.f90
View File

@ -24,9 +24,10 @@ SUBROUTINE dvqhub_barepsi_us2 (ik, dvqhbar, dvqhbar_orth, dvqhbar_orth_lm)
USE ions_base, ONLY : nat, ityp, ntyp => nsp
USE klist, ONLY : xk, ngk, igk_k
USE ldaU, ONLY : U_projection, Hubbard_l, is_hubbard, Hubbard_J0, offsetU, nwfcU
USE ldaU_ph, ONLY : wfcatomk, wfcatomkpq, swfcatomk, swfcatomkpq, dwfcatomkpq, &
USE ldaU_ph, ONLY : wfcatomk, wfcatomkpq, dwfcatomkpq, &
sdwfcatomk, sdwfcatomkpq, dvkb, vkbkpq, dvkbkpq, &
proj1, proj2, effU
proj1, proj2
USE ldaU_lr, ONLY : effU, swfcatomk, swfcatomkpq
USE wvfct, ONLY : npwx, nbnd
USE uspp, ONLY : vkb, nkb, okvan, ofsbeta
USE qpoint, ONLY : nksq, ikks, ikqs

4
PHonon/PH/dynmat_hub_bare.f90
View File

@ -33,8 +33,8 @@ SUBROUTINE dynmat_hub_bare
USE ldaU, ONLY : Hubbard_l, is_hubbard, Hubbard_J0, &
Hubbard_lmax, offsetU, nwfcU
USE ldaU_ph, ONLY : dnsbare, wfcatomk, wfcatomkpq, dvkb, vkbkpq, &
dvkbkpq, proj1, projpb, projpdb, swfcatomk, &
effU, read_dns_bare, d2ns_type
dvkbkpq, proj1, projpb, projpdb, read_dns_bare, d2ns_type
USE ldaU_lr, ONLY : effU, swfcatomk
USE wavefunctions, ONLY : evc
USE units_lr, ONLY : iuwfc, lrwfc, iuatwfc, iuatswfc
USE uspp, ONLY : vkb, nkb, ofsbeta

7
PHonon/PH/dynmat_hub_scf.f90
View File

@ -43,8 +43,8 @@ SUBROUTINE dynmat_hub_scf (irr, nu_i0, nper)
USE kinds, ONLY : DP
USE ions_base, ONLY : nat, ityp, ntyp => nsp
USE ldaU, ONLY : Hubbard_l, is_hubbard, Hubbard_J0
USE ldaU_ph, ONLY : dnsbare, dnsbare_all_modes, dnsscf, &
dnsorth_cart, effU
USE ldaU_lr, ONLY : dnsscf, effU
USE ldaU_ph, ONLY : dnsbare, dnsbare_all_modes, dnsorth_cart
USE lsda_mod, ONLY : lsda, current_spin, isk, nspin
USE modes, ONLY : u, nmodes
USE dynmat, ONLY : dyn, dyn_rec, dyn_hub_scf
@ -53,8 +53,7 @@ SUBROUTINE dynmat_hub_scf (irr, nu_i0, nper)
USE wvfct, ONLY : npwx, nbnd
USE control_lr, ONLY : lgamma
USE control_ph, ONLY : rec_code_read
USE units_ph, ONLY : iudwf, lrdwf
USE units_lr, ONLY : iuwfc, lrwfc
USE units_lr, ONLY : iuwfc, lrwfc, iudwf, lrdwf
USE wavefunctions, ONLY : evc
USE klist, ONLY : wk, lgauss, ltetra, ngk, igk_k
USE uspp, ONLY : okvan

6
PHonon/PH/ef_shift.f90
View File

@ -37,8 +37,7 @@ subroutine ef_shift (drhoscf, ldos, ldoss, dos_ef, irr, npe, flag)
USE qpoint, ONLY : nksq
USE control_lr, ONLY : nbnd_occ
USE control_ph, ONLY : lgamma_gamma
USE units_ph, ONLY : lrdwf, iudwf
USE units_lr, ONLY : iuwfc, lrwfc
USE units_lr, ONLY : iuwfc, lrwfc, lrdwf, iudwf
USE eqv, ONLY : dpsi
USE modes, ONLY : npert
USE mp_bands, ONLY : intra_bgrp_comm
@ -191,8 +190,7 @@ subroutine ef_shift_paw (drhoscf, dbecsum, ldos, ldoss, becsum1, &
USE control_lr, ONLY : nbnd_occ
USE control_ph, ONLY : lgamma_gamma
USE noncollin_module, ONLY : noncolin, npol, nspin_lsda, nspin_mag
USE units_ph, ONLY : lrdwf, iudwf
USE units_lr, ONLY : iuwfc, lrwfc
USE units_lr, ONLY : iuwfc, lrwfc, iudwf, lrdwf
USE eqv, ONLY : dpsi
USE modes, ONLY : npert
USE mp_bands, ONLY : intra_bgrp_comm

3
PHonon/PH/el_opt.f90
View File

@ -21,7 +21,8 @@ subroutine el_opt
USE symme, ONLY : symmatrix3
USE qpoint, ONLY : nksq
USE wvfct, ONLY : nbnd, npwx
USE units_ph, ONLY : iudrho, lrdrho, lrdwf, iudwf
USE units_ph, ONLY : iudrho, lrdrho
USE units_lr, ONLY : lrdwf, iudwf
USE control_lr, ONLY : nbnd_occ
USE buffers, ONLY : get_buffer
USE ph_restart, ONLY : ph_writefile

3
PHonon/PH/elphon.f90
View File

@ -352,7 +352,8 @@ SUBROUTINE elphel (irr, npe, imode0, dvscfins)
USE control_lr, ONLY : lgamma
USE fft_helper_subroutines
USE ldaU, ONLY : lda_plus_u, Hubbard_lmax
USE ldaU_ph, ONLY : dnsscf_all_modes, dnsscf
USE ldaU_lr, ONLY : dnsscf
USE ldaU_ph, ONLY : dnsscf_all_modes
USE io_global, ONLY : ionode, ionode_id
USE io_files, ONLY : seqopn
USE lrus, ONLY : becp1, int3_nc

25
PHonon/PH/make.depend
View File

@ -48,24 +48,6 @@ addcore.o : ../../Modules/kind.o
addcore.o : ../../Modules/recvec.o
addcore.o : ../../upflib/uspp.o
addcore.o : phcom.o
adddvhubscf.o : ../../LR_Modules/lrcom.o
adddvhubscf.o : ../../Modules/io_files.o
adddvhubscf.o : ../../Modules/ions_base.o
adddvhubscf.o : ../../Modules/kind.o
adddvhubscf.o : ../../Modules/mp_bands.o
adddvhubscf.o : ../../Modules/wavefunctions.o
adddvhubscf.o : ../../PW/src/buffers.o
adddvhubscf.o : ../../PW/src/ldaU.o
adddvhubscf.o : ../../PW/src/pwcom.o
adddvhubscf.o : ../../UtilXlib/mp.o
adddvhubscf.o : phcom.o
adddvscf_ph_mag.o : ../../LR_Modules/lrcom.o
adddvscf_ph_mag.o : ../../Modules/becmod.o
adddvscf_ph_mag.o : ../../Modules/ions_base.o
adddvscf_ph_mag.o : ../../Modules/kind.o
adddvscf_ph_mag.o : ../../Modules/noncol.o
adddvscf_ph_mag.o : ../../PW/src/pwcom.o
adddvscf_ph_mag.o : ../../upflib/uspp.o
addnlcc.o : ../../LR_Modules/lrcom.o
addnlcc.o : ../../Modules/cell_base.o
addnlcc.o : ../../Modules/fft_base.o
@ -252,7 +234,6 @@ chi_test.o : ../../Modules/kind.o
chi_test.o : ../../Modules/wavefunctions.o
chi_test.o : ../../PW/src/buffers.o
chi_test.o : ../../PW/src/pwcom.o
chi_test.o : phcom.o
chi_test.o : ramanm.o
clean_pw_ph.o : ../../LR_Modules/lrcom.o
clean_pw_ph.o : ../../Modules/kind.o
@ -437,7 +418,6 @@ dhdrhopsi.o : ../../PW/src/init_us_2.o
dhdrhopsi.o : ../../PW/src/pwcom.o
dhdrhopsi.o : ../../UtilXlib/mp.o
dhdrhopsi.o : ../../upflib/uspp.o
dhdrhopsi.o : phcom.o
dhdrhopsi.o : ramanm.o
dielec.o : ../../LR_Modules/lrcom.o
dielec.o : ../../Modules/cell_base.o
@ -1630,12 +1610,11 @@ set_irr_sym.o : phcom.o
setlocq.o : ../../Modules/constants.o
setlocq.o : ../../Modules/kind.o
setlocq.o : ../../PW/src/Coul_cut_2D.o
solve_e.o : ../../FFTXlib/fft_helper_subroutines.o
solve_e.o : ../../FFTXlib/fft_interfaces.o
solve_e.o : ../../LR_Modules/apply_dpot_mod.o
solve_e.o : ../../LR_Modules/dv_of_drho.o
solve_e.o : ../../LR_Modules/lrcom.o
solve_e.o : ../../Modules/cell_base.o
solve_e.o : ../../LR_Modules/response_kernels.o
solve_e.o : ../../Modules/check_stop.o
solve_e.o : ../../Modules/fft_base.o
solve_e.o : ../../Modules/io_files.o
@ -1711,11 +1690,11 @@ solve_e_nscf.o : ../../Modules/wavefunctions.o
solve_e_nscf.o : ../../PW/src/buffers.o
solve_e_nscf.o : ../../PW/src/pwcom.o
solve_e_nscf.o : phcom.o
solve_linter.o : ../../FFTXlib/fft_helper_subroutines.o
solve_linter.o : ../../FFTXlib/fft_interfaces.o
solve_linter.o : ../../LR_Modules/apply_dpot_mod.o
solve_linter.o : ../../LR_Modules/dv_of_drho.o
solve_linter.o : ../../LR_Modules/lrcom.o
solve_linter.o : ../../LR_Modules/response_kernels.o
solve_linter.o : ../../Modules/cell_base.o
solve_linter.o : ../../Modules/check_stop.o
solve_linter.o : ../../Modules/fft_base.o

6
PHonon/PH/openfilq.f90
View File

@ -14,9 +14,9 @@ SUBROUTINE openfilq()
!
USE kinds, ONLY : DP
USE control_flags, ONLY : io_level, modenum
USE units_ph, ONLY : iudwf, iubar, iucom, iudvkb3, &
USE units_ph, ONLY : iubar, iucom, iudvkb3, &
iudrhous, iuebar, iudrho, iudyn, iudvscf, &
lrdwf, lrbar, lrcom, lrdvkb3, &
lrbar, lrcom, lrdvkb3, &
lrdrhous, lrebar, lrdrho, lint3paw, iuint3paw, &
iundnsscf, iudvpsi, lrdvpsi, iugauge
USE units_lr, ONLY : iuwfc, lrwfc
@ -47,7 +47,7 @@ SUBROUTINE openfilq()
USE dfile_autoname, ONLY : dfile_name
USE qpoint, ONLY : xq
USE control_lr, ONLY : lgamma
USE units_lr, ONLY : iuatwfc, iuatswfc
USE units_lr, ONLY : iuatwfc, iuatswfc, iudwf, lrdwf
USE modes, ONLY : nmodes
USE ldaU, ONLY : lda_plus_u, Hubbard_lmax, nwfcU
USE ldaU_ph, ONLY : dnsscf_all_modes

2
PHonon/PH/ph_restart.f90
View File

@ -315,7 +315,7 @@ MODULE ph_restart
END DO
CALL xmlw_closetag( )
ENDDO
CALL xmlw_closetag( )
! Note: Root tag closed by routine ph_writefile
RETURN
END SUBROUTINE write_el_phon

7
PHonon/PH/phcom.f90
View File

@ -303,9 +303,7 @@ MODULE units_ph
lrbar, & ! length of the DV_{bare}
iuebar, & ! unit with the part DV_{bare} for the electric field
lrebar, & ! length of the DV_{bare} fro the electric field
iudwf, & ! unit with D psi
iupsir, & ! unit with evc in real space
lrdwf, & ! length of D psi record
iudrhous, lrdrhous, &
iudyn, & ! the unit for the dynamical matrix
iupdyn, & ! the unit for the partial dynamical matrix
@ -390,12 +388,10 @@ MODULE ldaU_ph
!
! atomic wfc's at k
COMPLEX(DP), ALLOCATABLE, TARGET :: wfcatomk(:,:), & ! atomic wfc at k
swfcatomk(:,:), & ! S * atomic wfc at k
dwfcatomk(:,:,:), & ! derivative of atomic wfc at k
sdwfcatomk(:,:) ! S * derivative of atomic wfc at k
! atomic wfc's at k+q
COMPLEX(DP), POINTER :: wfcatomkpq(:,:), & ! atomic wfc at k+q
swfcatomkpq(:,:), & ! S * atomic wfc at k+q
dwfcatomkpq(:,:,:), & ! derivative of atomic wfc at k+q
sdwfcatomkpq(:,:) ! S * derivative of atomic wfc at k+q
!
@ -406,7 +402,6 @@ MODULE ldaU_ph
! Various arrays for the response occupation matrix
COMPLEX(DP), ALLOCATABLE :: dnsbare(:,:,:,:,:,:), & ! bare derivative of ns
dnsbare_all_modes(:,:,:,:,:), & ! bare derivative of ns for all modes
dnsscf(:,:,:,:,:), & ! SCF derivative of ns
dnsscf_all_modes(:,:,:,:,:), & ! SCF derivative of ns for all modes
dnsorth(:,:,:,:,:), & ! valence component of dns
dnsorth_cart(:,:,:,:,:,:) ! same as above, but in cart. coordinates
@ -420,8 +415,6 @@ MODULE ldaU_ph
! projpdb = <psi|dbeta>
!
!
REAL(DP) :: effU(ntypx)
! effective Hubbard parameter: effU = Hubbard_U - Hubbard_J0
LOGICAL :: read_dns_bare
! if .true. read the first bare derivative of ns from file
CHARACTER(LEN=4) :: d2ns_type

3
PHonon/PH/phescf.f90
View File

@ -28,7 +28,8 @@ SUBROUTINE phescf()
USE lrus, ONLY : int3, int3_nc, int3_paw
USE control_lr, ONLY : lrpa
USE ldaU, ONLY : lda_plus_u, Hubbard_lmax
USE ldaU_ph, ONLY : dnsscf, dnsscf_all_modes
USE ldaU_lr, ONLY : dnsscf
USE ldaU_ph, ONLY : dnsscf_all_modes
USE control_flags, ONLY : iverbosity
USE write_hub
!

2
PHonon/PH/phq_setup.f90
View File

@ -99,7 +99,7 @@ subroutine phq_setup
USE nc_mag_aux, ONLY : deeq_nc_save
USE control_lr, ONLY : lgamma
USE ldaU, ONLY : lda_plus_u, Hubbard_U, Hubbard_J0
USE ldaU_ph, ONLY : effU
USE ldaU_lr, ONLY : effU
USE constants, ONLY : rytoev
USE dvscf_interpolate, ONLY : ldvscf_interpolate, dvscf_interpol_setup
USE ahc, ONLY : elph_ahc, elph_ahc_setup

2
PHonon/PH/phq_summary.f90
View File

@ -44,7 +44,7 @@ subroutine phq_summary
USE qpoint, ONLY : xq
USE lr_symm_base, ONLY : irotmq, minus_q, nsymq
USE constants, ONLY : rytoev
USE ldaU_ph, ONLY : effU
USE ldaU_lr, ONLY : effU
USE ldaU, ONLY : lda_plus_u, Hubbard_U, Hubbard_J0, &
lda_plus_u_kind, is_hubbard

3
PHonon/PH/phqscf.f90
View File

@ -36,7 +36,8 @@ SUBROUTINE phqscf
USE eqv, ONLY : drhoscfs
USE dynmat, ONLY : dyn_hub_scf
USE ldaU, ONLY : lda_plus_u, Hubbard_lmax
USE ldaU_ph, ONLY : dnsscf, dnsscf_all_modes
USE ldaU_lr, ONLY : dnsscf
USE ldaU_ph, ONLY : dnsscf_all_modes
USE units_ph, ONLY : iundnsscf
USE control_flags, ONLY : iverbosity
USE write_hub

3
PHonon/PH/polariz.f90
View File

@ -21,7 +21,8 @@ subroutine polariz ( iw, iu )
USE wvfct, ONLY : npwx
USE kinds, ONLY : DP
USE control_lr, ONLY : nbnd_occ
USE units_ph, ONLY : lrdwf, iudwf, lrebar, iuebar
USE units_ph, ONLY : lrebar, iuebar
USE units_lr, ONLY : iudwf, lrdwf
USE buffers, ONLY : get_buffer
USE freq_ph, ONLY : polar, done_iu, comp_iu
USE eqv, ONLY : dpsi, dvpsi

3
PHonon/PH/raman_mat.f90
View File

@ -26,8 +26,7 @@ subroutine raman_mat
USE wvfct, ONLY : npwx, nbnd
USE wavefunctions, ONLY: evc
USE phus, ONLY : alphap
USE units_ph, ONLY : lrdwf, iudwf
USE units_lr, ONLY : iuwfc, lrwfc
USE units_lr, ONLY : iuwfc, lrwfc, iudwf, lrdwf
USE ramanm, ONLY : ramtns, jab, a1j, a2j, lrd2w, iud2w
USE lrus, ONLY : becp1

265
PHonon/PH/solve_e.f90
View File

@ -22,31 +22,33 @@ subroutine solve_e
! e) computes Delta rho, Delta V_{SCF} and symmetrizes them
! f) If lda_plus_u=.true. compute also the response occupation
! matrices dnsscf
! Step b, c, d are done in sternheimer_kernel.
!
USE kinds, ONLY : DP
USE ions_base, ONLY : nat
USE io_global, ONLY : stdout, ionode
USE io_files, ONLY : prefix, diropn
USE cell_base, ONLY : tpiba2
USE klist, ONLY : ltetra, lgauss, xk, wk, ngk, igk_k
USE gvect, ONLY : g
USE io_files, ONLY : diropn
USE mp, ONLY : mp_sum
USE mp_pools, ONLY : inter_pool_comm
USE mp_bands, ONLY : intra_bgrp_comm
USE klist, ONLY : ltetra, lgauss, xk, ngk, igk_k
USE gvecs, ONLY : doublegrid
USE fft_base, ONLY : dfftp, dffts
USE lsda_mod, ONLY : lsda, nspin, current_spin, isk
USE lsda_mod, ONLY : nspin
USE spin_orb, ONLY : domag
USE wvfct, ONLY : nbnd, npwx, et
USE wvfct, ONLY : nbnd, npwx
USE check_stop, ONLY : check_stop_now
USE buffers, ONLY : get_buffer, save_buffer
USE buffers, ONLY : get_buffer
USE wavefunctions, ONLY : evc
USE uspp, ONLY : okvan, vkb
USE uspp_param, ONLY : upf, nhm
USE uspp_param, ONLY : nhm
USE noncollin_module, ONLY : noncolin, npol, nspin_mag
USE scf, ONLY : rho
USE paw_variables, ONLY : okpaw
USE paw_onecenter, ONLY : paw_dpotential
USE paw_symmetry, ONLY : paw_desymmetrize
USE units_ph, ONLY : lrdwf, iudwf, lrdrho, iudrho
USE units_ph, ONLY : lrdrho, iudrho, lrebar, iuebar
USE units_lr, ONLY : iuwfc, lrwfc
USE output, ONLY : fildrho
USE control_ph, ONLY : ext_recover, rec_code, &
@ -54,57 +56,46 @@ subroutine solve_e
alpha_mix, lgamma_gamma, niter_ph, &
flmixdpot, rec_code_read
USE recover_mod, ONLY : read_rec, write_rec
USE mp_pools, ONLY : inter_pool_comm
USE mp_bands, ONLY : intra_bgrp_comm, ntask_groups
USE mp, ONLY : mp_sum
USE lrus, ONLY : int3_paw
USE qpoint, ONLY : nksq, ikks, ikqs
USE eqv, ONLY : dpsi, dvpsi
USE control_lr, ONLY : nbnd_occ, lgamma
USE dv_of_drho_lr
USE fft_helper_subroutines
USE qpoint, ONLY : nksq, ikks
USE control_lr, ONLY : lgamma
USE dv_of_drho_lr, ONLY : dv_of_drho
USE fft_interfaces, ONLY : fft_interpolate
USE ldaU, ONLY : lda_plus_u
USE apply_dpot_mod, ONLY : apply_dpot_allocate, apply_dpot_deallocate, &
apply_dpot_bands
USE apply_dpot_mod, ONLY : apply_dpot_allocate, apply_dpot_deallocate
USE response_kernels, ONLY : sternheimer_kernel
USE uspp_init, ONLY : init_us_2
!
IMPLICIT NONE
!
LOGICAL :: exst
!!
LOGICAL :: all_conv
!! True if sternheimer_kernel is converged at all k points and perturbations
INTEGER :: ikk, npw, kter, iter0, ipol, iter, ik, is, ndim
!! counters
REAL(DP) :: thresh
!! convergence threshold
REAL(DP) :: averlt
!! average number of iterations
REAL(DP) :: dr2
!! self-consistency error
REAL(DP) :: tcpu, get_clock
!! timing variables
implicit none
real(DP) :: thresh, anorm, averlt, dr2
! thresh: convergence threshold
! anorm : the norm of the error
! averlt: average number of iterations
! dr2 : self-consistency error
real(DP), allocatable :: h_diag (:,:)
! h_diag: diagonal part of the Hamiltonian
complex(DP) , allocatable, target :: &
dvscfin (:,:,:) ! change of the scf potential (input)
complex(DP) , pointer :: &
dvscfins (:,:,:) ! change of the scf potential (smooth)
complex(DP) , allocatable :: &
dvscfout (:,:,:), & ! change of the scf potential (output)
dbecsum(:,:,:,:), & ! the becsum with dpsi
dbecsum_nc(:,:,:,:,:), & ! the becsum with dpsi
mixin(:), mixout(:), & ! auxiliary for paw mixing
ps (:,:), &
aux2(:,:)
logical :: conv_root, exst
! conv_root: true if linear system is converged
integer :: npw, npwq
integer :: kter, iter0, ipol, ibnd, iter, lter, ik, ig, is, nrec, ndim, ios
! counters
integer :: ltaver, lintercall
integer :: ikk, ikq
real(DP) :: tcpu, get_clock
! timing variables
external ch_psi_all, cg_psi
COMPLEX(DP), ALLOCATABLE, TARGET :: dvscfin (:,:,:)
!! change of the scf potential (input)
COMPLEX(DP), POINTER :: dvscfins (:,:,:)
!! change of the scf potential (smooth)
COMPLEX(DP), ALLOCATABLE :: dvscfout (:,:,:)
!! change of the scf potential (output)
COMPLEX(DP), ALLOCATABLE :: dbecsum(:,:,:,:)
!! the becsum with dpsi
COMPLEX(DP), ALLOCATABLE :: dbecsum_nc(:,:,:,:,:)
!! the becsum with dpsi
COMPLEX(DP), ALLOCATABLE :: mixin(:), mixout(:)
!! auxiliary for paw mixing
!
call start_clock ('solve_e')
!
! This routine is task group aware
@ -120,12 +111,10 @@ subroutine solve_e
IF (okpaw) THEN
ALLOCATE (mixin(dfftp%nnr*nspin_mag*3+(nhm*(nhm+1)*nat*nspin_mag*3)/2) )
ALLOCATE (mixout(dfftp%nnr*nspin_mag*3+(nhm*(nhm+1)*nat*nspin_mag*3)/2) )
mixin=(0.0_DP,0.0_DP)
ENDIF
allocate (dbecsum( nhm*(nhm+1)/2, nat, nspin_mag, 3))
IF (noncolin) allocate (dbecsum_nc (nhm, nhm, nat, nspin, 3))
allocate (h_diag(npwx*npol, nbnd))
allocate (aux2(npwx*npol, nbnd))
IF (okpaw) mixin=(0.0_DP,0.0_DP)
CALL apply_dpot_allocate()
if (rec_code_read == -20.AND.ext_recover) then
@ -159,135 +148,57 @@ subroutine solve_e
if ( (lgauss .or. ltetra) .or..not.lgamma) call errore ('solve_e', &
'called in the wrong case', 1)
!
! Compute P_c^+ x psi for all polarization and k points and store in buffer
!
DO ik = 1, nksq
DO ipol = 1, 3
ikk = ikks(ik)
npw = ngk(ikk)
!
! reads unperturbed wavefunctions psi_k in G_space, for all bands
!
IF (nksq > 1) THEN
CALL get_buffer(evc, lrwfc, iuwfc, ikk)
ENDIF
!
CALL init_us_2(npw, igk_k(1, ikk), xk(1, ikk), vkb)
!
! computes P_c^+ x psi_kpoint, written to buffer iuebar.
!
CALL dvpsi_e(ik, ipol)
!
ENDDO ! ipol
ENDDO ! ik
!
! The outside loop is over the iterations
!
do kter = 1, niter_ph
!
FLUSH( stdout )
iter = kter + iter0
ltaver = 0
lintercall = 0
!
dvscfout(:,:,:)=(0.d0,0.d0)
dbecsum(:,:,:,:)=(0.d0,0.d0)
IF (noncolin) dbecsum_nc=(0.d0,0.d0)
dvscfout = (0.d0,0.d0)
dbecsum = (0.d0,0.d0)
IF (noncolin) dbecsum_nc = (0.d0,0.d0)
!
! DFPT+U: at each iteration calculate dnsscf,
! i.e. the scf variation of the occupation matrix ns.
!
IF (lda_plus_u .AND. (iter.NE.1)) &
CALL dnsq_scf (3, .false., 0, 1, .false.)
IF (lda_plus_u .AND. (iter /= 1)) CALL dnsq_scf(3, .false., 0, 1, .false.)
!
do ik = 1, nksq
ikk=ikks(ik)
ikq=ikqs(ik)
!
npw = ngk(ikk)
npwq= npw ! q=0 always in this routine
!
if (lsda) current_spin = isk (ikk)
!
! reads unperturbed wavefunctions psi_k in G_space, for all bands
!
if (nksq.ge.1) THEN
call get_buffer (evc, lrwfc, iuwfc, ikk)
end if
!
! compute beta functions and kinetic energy for k-point ik
! needed by h_psi, called by ch_psi_all, called by cgsolve_all
!
CALL init_us_2 (npw, igk_k(1,ikk), xk (1, ikk), vkb)
CALL g2_kin(ikk)
!
! compute preconditioning matrix h_diag used by cgsolve_all
!
CALL h_prec (ik, evc, h_diag)
!
do ipol = 1, 3
!
! computes/reads P_c^+ x psi_kpoint into dvpsi array
!
call dvpsi_e (ik, ipol)
!
if (iter > 1) then
!
! calculates dvscf_q*psi_k in G_space, for all bands, k=kpoint
! dvscf_q from previous iteration (mix_potential)
!
CALL apply_dpot_bands(ik, nbnd_occ(ikk), dvscfins(:, :, ipol), evc, aux2)
dvpsi = dvpsi + aux2
!
! In the case of US pseudopotentials there is an additional
! selfconsist term which comes from the dependence of D on
! V_{eff} on the bare change of the potential
!
call adddvscf(ipol,ik)
!
! DFPT+U: add to dvpsi the scf part of the response
! Hubbard potential dV_hub
!
if (lda_plus_u) call adddvhubscf (ipol, ik)
!
endif
!
! Orthogonalize dvpsi to valence states: ps = <evc|dvpsi>
!
CALL orthogonalize(dvpsi, evc, ikk, ikk, dpsi, npwq, .false.)
!
if (iter == 1) then
!
! At the first iteration dpsi and dvscfin are set to zero,
!
dpsi(:,:)=(0.d0,0.d0)
dvscfin(:,:,:)=(0.d0,0.d0)
!
! starting threshold for the iterative solution of the linear
! system
!
thresh = 1.d-2
if (lnoloc) thresh = 1.d-5
else
! starting value for delta_psi is read from iudwf
!
nrec = (ipol - 1) * nksq + ik
call get_buffer (dpsi, lrdwf, iudwf, nrec)
!
! threshold for iterative solution of the linear system
!
thresh = min (0.1d0 * sqrt (dr2), 1.0d-2)
endif
!
! iterative solution of the linear system (H-e)*dpsi=dvpsi
! dvpsi=-P_c+ (dvbare+dvscf)*psi , dvscf fixed.
!
conv_root = .true.
call cgsolve_all (ch_psi_all,cg_psi,et(1,ikk),dvpsi,dpsi, &
h_diag,npwx,npw,thresh,ik,lter,conv_root,anorm,nbnd_occ(ikk),npol)
ltaver = ltaver + lter
lintercall = lintercall + 1
if (.not.conv_root) WRITE( stdout, "(5x,'kpoint',i4,' ibnd',i4, &
& ' solve_e: root not converged ',es10.3)") ik &
&, ibnd, anorm
!
! writes delta_psi on iunit iudwf, k=kpoint,
!
nrec = (ipol - 1) * nksq + ik
call save_buffer(dpsi, lrdwf, iudwf, nrec)
!
! calculates dvscf, sum over k => dvscf_q_ipert
!
IF (noncolin) THEN
call incdrhoscf_nc(dvscfout(1,1,ipol),wk(ikk),ik, &
dbecsum_nc(1,1,1,1,ipol), dpsi, 1.0d0)
ELSE
call incdrhoscf (dvscfout(1,current_spin,ipol), wk(ikk), &
ik, dbecsum(1,1,current_spin,ipol), dpsi)
ENDIF
enddo ! on polarizations
enddo ! on k points
! set threshold for the iterative solution of the linear system
!
IF (iter == 1) THEN
thresh = 1.d-2
IF (lnoloc) thresh = 1.d-5
ELSE
thresh = MIN(0.1d0 * SQRT(dr2), 1.0d-2)
ENDIF
!
! Compute dvscfout, the charge density response to the total potential
!
CALL sternheimer_kernel(iter==1, .FALSE., 3, lrebar, iuebar, thresh, dvscfins, &
all_conv, averlt, dvscfout, dbecsum, dbecsum_nc)
!
! The calculation of dbecsum is distributed across processors
! (see addusdbec) - we sum over processors the contributions
@ -374,8 +285,6 @@ subroutine solve_e
call newdq(dvscfin,3)
averlt = DBLE (ltaver) / DBLE (lintercall)
tcpu = get_clock ('PHONON')
WRITE( stdout, '(/,5x," iter # ",i3," total cpu time :",f8.1, &
& " secs av.it.: ",f5.1)') iter, tcpu, averlt
@ -404,7 +313,6 @@ subroutine solve_e
155 continue
!
CALL apply_dpot_deallocate()
deallocate (h_diag)
deallocate (dbecsum)
deallocate (dvscfout)
IF (okpaw) THEN
@ -414,7 +322,6 @@ subroutine solve_e
if (doublegrid) deallocate (dvscfins)
deallocate (dvscfin)
if (noncolin) deallocate(dbecsum_nc)
deallocate(aux2)
call stop_clock ('solve_e')
return

6
PHonon/PH/solve_e_fpol.f90
View File

@ -16,7 +16,7 @@ subroutine solve_e_fpol ( iw )
! a) computes the bare potential term x | psi >
! b) adds to it the screening term Delta V_{SCF} | psi >
! c) applies P_c^+ (orthogonalization to valence states)
! d) calls cgsolve_all to solve the linear system
! d) calls gmressolve_all to solve the linear system
! e) computes Delta rho, Delta V_{SCF} and symmetrizes them
!
USE kinds, ONLY : DP
@ -42,8 +42,8 @@ subroutine solve_e_fpol ( iw )
rec_code, flmixdpot
USE output, ONLY : fildrho
USE qpoint, ONLY : nksq
USE units_ph, ONLY : lrdwf, iudwf, iudrho, lrdrho
USE units_lr, ONLY : iuwfc, lrwfc
USE units_ph, ONLY : iudrho, lrdrho
USE units_lr, ONLY : iuwfc, lrwfc, iudwf, lrdwf
USE mp_pools, ONLY : inter_pool_comm
USE mp_bands, ONLY : intra_bgrp_comm
USE mp, ONLY : mp_sum

3
PHonon/PH/solve_e_nscf.f90
View File

@ -25,8 +25,9 @@ subroutine solve_e_nscf( avg_iter, thresh, ik, ipol, dvscfs, auxr )
USE wvfct, ONLY : et
USE wavefunctions, ONLY : evc
USE eqv, ONLY : dpsi, dvpsi
USE units_ph, ONLY : this_pcxpsi_is_on_file, lrdwf, iudwf
USE units_ph, ONLY : this_pcxpsi_is_on_file
USE qpoint, ONLY : nksq
USE units_lr, ONLY : lrdwf, iudwf
USE control_lr, ONLY : nbnd_occ
implicit none

381
PHonon/PH/solve_linter.f90
View File

@ -25,29 +25,29 @@ SUBROUTINE solve_linter (irr, imode0, npe, drhoscf)
! e) computes Delta rho, Delta V_{SCF} and symmetrizes them
! f) If lda_plus_u=.true. compute also the response occupation
! matrices dnsscf
! g) (Introduced in February 2020) If noncolin=.true. and domag=.true.
! the linear system is solved twice (nsolv = 2, the case
! isolv = 2 needs the time-reversed wave functions). For the
! theoretical background, please refer to Phys. Rev. B 100,
! g) (Introduced in February 2020) If noncolin=.true. and domag=.true.
! the linear system is solved twice (nsolv = 2, the case
! isolv = 2 needs the time-reversed wave functions). For the
! theoretical background, please refer to Phys. Rev. B 100,
! 045115 (2019)
! Step b, c, d are done inside sternheimer_kernel.
!
USE kinds, ONLY : DP
USE ions_base, ONLY : nat, ntyp => nsp, ityp
USE ions_base, ONLY : nat
USE io_global, ONLY : stdout, ionode
USE io_files, ONLY : prefix, diropn
USE check_stop, ONLY : check_stop_now
USE wavefunctions, ONLY : evc
USE cell_base, ONLY : at
USE klist, ONLY : ltetra, lgauss, xk, wk, ngk, igk_k
USE gvect, ONLY : g
USE klist, ONLY : ltetra, lgauss, xk, ngk, igk_k
USE gvecs, ONLY : doublegrid
USE fft_base, ONLY : dfftp, dffts
USE lsda_mod, ONLY : lsda, nspin, current_spin, isk
USE spin_orb, ONLY : domag
USE wvfct, ONLY : nbnd, npwx, et
USE wvfct, ONLY : nbnd, npwx
USE scf, ONLY : rho, vrs
USE uspp, ONLY : okvan, vkb, deeq_nc
USE uspp_param, ONLY : upf, nhm
USE uspp_param, ONLY : nhm
USE noncollin_module, ONLY : noncolin, npol, nspin_mag
USE paw_variables, ONLY : okpaw
USE paw_onecenter, ONLY : paw_dpotential
@ -59,12 +59,12 @@ SUBROUTINE solve_linter (irr, imode0, npe, drhoscf)
where_rec, flmixdpot, ext_recover
USE el_phon, ONLY : elph
USE uspp, ONLY : nlcc_any
USE units_ph, ONLY : iudrho, lrdrho, iudwf, lrdwf, iubar, lrbar, &
USE units_ph, ONLY : iudrho, lrdrho, iubar, lrbar, &
iudvscf, iuint3paw, lint3paw
USE units_lr, ONLY : iuwfc, lrwfc
USE output, ONLY : fildrho, fildvscf
USE phus, ONLY : becsumort, alphap, int1_nc
USE modes, ONLY : npertx, npert, u, t, tmq
USE modes, ONLY : npertx, u, t, tmq
USE recover_mod, ONLY : read_rec, write_rec
! used to write fildrho:
USE dfile_autoname, ONLY : dfile_name
@ -76,17 +76,16 @@ SUBROUTINE solve_linter (irr, imode0, npe, drhoscf)
USE efermi_shift, ONLY : ef_shift, ef_shift_paw, def
USE lrus, ONLY : int3_paw, becp1, int3_nc
USE lr_symm_base, ONLY : irotmq, minus_q, nsymq, rtau
USE eqv, ONLY : dvpsi, dpsi, evq
USE eqv, ONLY : dvpsi
USE qpoint, ONLY : xq, nksq, ikks, ikqs
USE qpoint_aux, ONLY : ikmks, ikmkmqs, becpt, alphapt
USE control_lr, ONLY : nbnd_occ, lgamma
USE qpoint_aux, ONLY : ikmks, becpt, alphapt
USE control_lr, ONLY : lgamma
USE dv_of_drho_lr, ONLY : dv_of_drho
USE fft_helper_subroutines
USE fft_interfaces, ONLY : fft_interpolate
USE ldaU, ONLY : lda_plus_u
USE nc_mag_aux, ONLY : int1_nc_save, deeq_nc_save, int3_save
USE apply_dpot_mod, ONLY : apply_dpot_allocate, apply_dpot_deallocate, &
apply_dpot_bands
USE apply_dpot_mod, ONLY : apply_dpot_allocate, apply_dpot_deallocate
USE response_kernels, ONLY : sternheimer_kernel
USE uspp_init, ONLY : init_us_2
implicit none
@ -99,15 +98,11 @@ SUBROUTINE solve_linter (irr, imode0, npe, drhoscf)
complex(DP) :: drhoscf (dfftp%nnr, nspin_mag, npe)
! output: the change of the scf charge
real(DP) , allocatable :: h_diag (:,:)
! h_diag: diagonal part of the Hamiltonian
real(DP) :: thresh, anorm, averlt, dr2, rsign
real(DP) :: thresh, averlt, dr2
! thresh: convergence threshold
! anorm : the norm of the error
! averlt: average number of iterations
! dr2 : self-consistency error
! rsign : sign or the term in the magnetization
real(DP) :: dos_ef, weight, aux_avg (2)
real(DP) :: dos_ef
! Misc variables for metals
! dos_ef: density of states at Ef
@ -128,38 +123,30 @@ SUBROUTINE solve_linter (irr, imode0, npe, drhoscf)
! drhoc: response core charge density
REAL(DP), allocatable :: becsum1(:,:,:)
logical :: conv_root, & ! true if linear system is converged
exst, & ! used to open the recover file
lmetq0 ! true if xq=(0,0,0) in a metal
LOGICAL :: all_conv
!! True if sternheimer_kernel is converged at all k points and perturbations
logical :: exst, & ! used to open the recover file
lmetq0, & ! true if xq=(0,0,0) in a metal
first_iter ! true if first iteration where induced rho is not yet calculated
integer :: kter, & ! counter on iterations
iter0, & ! starting iteration
ipert, & ! counter on perturbations
ibnd, & ! counter on bands
iter, & ! counter on iterations
lter, & ! counter on iterations of linear system
ltaver, & ! average counter
lintercall, & ! average number of calls to cgsolve_all
ik, ikk, & ! counter on k points
ikq, & ! counter on k+q points
ig, & ! counter on G vectors
ndim, &
is, & ! counter on spin polarizations
nrec, & ! the record number for dvpsi and dpsi
ios, & ! integer variable for I/O control
ipol, & ! counter on polarization
mode, & ! mode index
mode, & ! mode index
isolv, & ! counter on linear systems
nsolv, & ! number of linear systems
ikmk, & ! index of mk
ikmkmq ! index of mk-mq
ikmk ! index of mk
integer :: npw, npwq
integer :: iq_dummy
real(DP) :: tcpu, get_clock ! timing variables
character(len=256) :: filename
external ch_psi_all, cg_psi
!
IF (rec_code_read > 20 ) RETURN
@ -188,7 +175,6 @@ SUBROUTINE solve_linter (irr, imode0, npe, drhoscf)
ALLOCATE(mixin(1))
ENDIF
IF (noncolin) allocate (dbecsum_nc (nhm,nhm, nat , nspin , npe, nsolv))
allocate (h_diag ( npwx*npol, nbnd))
allocate (aux2(npwx*npol, nbnd))
allocate (drhoc(dfftp%nnr))
IF (noncolin.AND.domag.AND.okvan) THEN
@ -245,220 +231,129 @@ SUBROUTINE solve_linter (irr, imode0, npe, drhoscf)
!
IF (iter0==-1000) iter0=0
!
! Compute dV_bare * psi and write to buffer iubar
!
DO ik = 1, nksq
!
ikk = ikks(ik)
ikq = ikqs(ik)
npw = ngk(ikk)
npwq = ngk(ikq)
!
IF (lsda) current_spin = isk(ikk)
!
! compute beta functions for k-point ikq
!
CALL init_us_2(npwq, igk_k(1, ikq), xk(1, ikq), vkb)
!
DO isolv = 1, nsolv
IF (isolv == 1) THEN
ikmk = ikks(ik)
ELSE
ikmk = ikmks(ik)
ENDIF
!
! read unperturbed wavefunctions psi(k) and psi(k+q)
!
IF (nksq > 1 .OR. nsolv == 2) THEN
CALL get_buffer(evc, lrwfc, iuwfc, ikmk)
ENDIF
!
DO ipert = 1, npe
mode = imode0 + ipert
nrec = (isolv-1) * npe * nksq + (ipert - 1) * nksq + ik
!
IF (isolv==1) THEN
CALL dvqpsi_us(ik, u(1, mode), .FAlSE., becp1, alphap)
!
! DFPT+U: At the first ph iteration the bare perturbed
! Hubbard potential dvbare_hub_q * psi_kpoint
! is calculated and added to dvpsi.
!
IF (lda_plus_u) CALL dvqhub_barepsi_us(ik, u(1, mode))
!
ELSE
IF (okvan) THEN
deeq_nc(:,:,:,:) = deeq_nc_save(:,:,:,:,2)
int1_nc(:,:,:,:,:) = int1_nc_save(:,:,:,:,:,2)
ENDIF
CALL dvqpsi_us(ik, u(1, mode), .FAlSE., becpt, alphapt)
IF (okvan) THEN
deeq_nc(:,:,:,:) = deeq_nc_save(:,:,:,:,1)
int1_nc(:,:,:,:,:) = int1_nc_save(:,:,:,:,:,1)
ENDIF
ENDIF
!
CALL save_buffer(dvpsi, lrbar, iubar, nrec)
!
ENDDO ! ipert
ENDDO ! isolv
ENDDO ! ik
!
! The outside loop is over the iterations
!
do kter = 1, niter_ph
!
iter = kter + iter0
ltaver = 0
lintercall = 0
!
drhoscf(:,:,:) = (0.d0, 0.d0)
dbecsum(:,:,:,:) = (0.d0, 0.d0)
first_iter = .NOT. (where_rec == 'solve_lint' .OR. iter > 1)
!
drhoscf = (0.d0, 0.d0)
dbecsum = (0.d0, 0.d0)
IF (noncolin) dbecsum_nc = (0.d0, 0.d0)
!
! DFPT+U: at each ph iteration calculate dnsscf,
! i.e. the scf variation of the occupation matrix ns.
!
IF (lda_plus_u .AND. (iter.NE.1)) &
CALL dnsq_scf (npe, lmetq0, imode0, irr, .true.)
IF (lda_plus_u .AND. (iter /= 1)) CALL dnsq_scf(npe, lmetq0, imode0, irr, .true.)
!
do ik = 1, nksq
! Start the loop on the two linear systems, one at B and one at -B
!
DO isolv = 1, nsolv
!
ikk = ikks(ik)
ikq = ikqs(ik)
npw = ngk(ikk)
npwq= ngk(ikq)
! change the sign of the magnetic field if required
!
if (lsda) current_spin = isk (ikk)
!
! compute beta functions and kinetic energy for k-point ikq
! needed by h_psi, called by ch_psi_all, called by cgsolve_all
!
CALL init_us_2 (npwq, igk_k(1,ikq), xk (1, ikq), vkb)
CALL g2_kin (ikq)
!
! Start the loop on the two linear systems, one at B and one at
! -B
!
DO isolv=1, nsolv
IF (isolv==2) THEN
ikmk = ikmks(ik)
ikmkmq = ikmkmqs(ik)
rsign=-1.0_DP
ELSE
ikmk=ikk
ikmkmq=ikq
rsign=1.0_DP
IF (isolv == 2) THEN
IF (.NOT. first_iter) THEN
dvscfins(:, 2:4, :) = -dvscfins(:, 2:4, :)
IF (okvan) int3_nc(:,:,:,:,:) = int3_save(:,:,:,:,:,2)
ENDIF
!
! read unperturbed wavefunctions psi(k) and psi(k+q)
!
if (nksq.gt.1.OR.nsolv==2) then
if (lgamma) then
call get_buffer (evc, lrwfc, iuwfc, ikmk)
else
call get_buffer (evc, lrwfc, iuwfc, ikmk)
call get_buffer (evq, lrwfc, iuwfc, ikmkmq)
end if
endif
!
! compute preconditioning matrix h_diag used by cgsolve_all
!
CALL h_prec (ik, evq, h_diag)
!
do ipert = 1, npe
mode = imode0 + ipert
nrec = (ipert - 1) * nksq + ik + (isolv-1) * npe * nksq
!
! and now adds the contribution of the self consistent term
!
if (where_rec =='solve_lint'.or.iter>1) then
!
! After the first iteration dvbare_q*psi_kpoint is read from file
!
call get_buffer (dvpsi, lrbar, iubar, nrec)
!
! calculates dvscf_q*psi_k in G_space, for all bands, k=kpoint
! dvscf_q from previous iteration (mix_potential)
!
call start_clock ('vpsifft')
!
! change the sign of the magnetic field if required
!
IF (isolv==2) THEN
dvscfins(:,2:4,ipert)=-dvscfins(:,2:4,ipert)
IF (okvan) int3_nc(:,:,:,:,ipert)=int3_save(:,:,:,:,ipert,2)
ENDIF
!
CALL apply_dpot_bands(ik, nbnd_occ(ikk), dvscfins(:, :, ipert), evc, aux2)
dvpsi = dvpsi + aux2
!
call stop_clock ('vpsifft')
!
! In the case of US pseudopotentials there is an additional
! selfconsist term which comes from the dependence of D on
! V_{eff} on the bare change of the potential
!
IF (isolv==1) THEN
call adddvscf_ph_mag (ipert, ik, becp1)
!
! DFPT+U: add to dvpsi the scf part of the response
! Hubbard potential dV_hub
!
if (lda_plus_u) call adddvhubscf (ipert, ik)
ELSE
call adddvscf_ph_mag (ipert, ik, becpt)
END IF
!
! reset the original magnetic field if it was changed
!
IF (isolv==2) THEN
dvscfins(:,2:4,ipert)=-dvscfins(:,2:4,ipert)
IF (okvan) int3_nc(:,:,:,:,ipert)=int3_save(:,:,:,:,ipert,1)
ENDIF
!
else
!
! At the first iteration dvbare_q*psi_kpoint is calculated
! and written to file.
!
IF (isolv==1) THEN
call dvqpsi_us (ik, u (1, mode),.false., becp1, alphap )
!
! DFPT+U: At the first ph iteration the bare perturbed
! Hubbard potential dvbare_hub_q * psi_kpoint
! is calculated and added to dvpsi.
!
if (lda_plus_u) call dvqhub_barepsi_us (ik, u(1,mode))
!
ELSE
IF (okvan) THEN
deeq_nc(:,:,:,:)=deeq_nc_save(:,:,:,:,2)
int1_nc(:,:,:,:,:)=int1_nc_save(:,:,:,:,:,2)
ENDIF
call dvqpsi_us (ik, u (1, mode),.false., becpt, alphapt)
IF (okvan) THEN
deeq_nc(:,:,:,:)=deeq_nc_save(:,:,:,:,1)
int1_nc(:,:,:,:,:)=int1_nc_save(:,:,:,:,:,1)
ENDIF
ENDIF
call save_buffer (dvpsi, lrbar, iubar, nrec)
!
endif
!
! Ortogonalize dvpsi to valence states: ps = <evq|dvpsi>
! Apply -P_c^+.
!
CALL orthogonalize(dvpsi, evq, ikmk, ikmkmq, dpsi, npwq, .false.)
!
if (where_rec=='solve_lint'.or.iter > 1) then
!
! starting value for delta_psi is read from iudwf
!
call get_buffer( dpsi, lrdwf, iudwf, nrec)
!
! threshold for iterative solution of the linear system
!
thresh = min (1.d-1 * sqrt (dr2), 1.d-2)
else
!
! At the first iteration dpsi and dvscfin are set to zero
!
dpsi(:,:) = (0.d0, 0.d0)
dvscfin (:, :, ipert) = (0.d0, 0.d0)
!
! starting threshold for iterative solution of the linear system
!
thresh = 1.0d-2
endif
!
! iterative solution of the linear system (H-eS)*dpsi=dvpsi,
! dvpsi=-P_c^+ (dvbare+dvscf)*psi , dvscf fixed.
!
IF (isolv==2) THEN
vrs(:,2:4)=-vrs(:,2:4)
IF (okvan) deeq_nc(:,:,:,:)=deeq_nc_save(:,:,:,:,2)
ENDIF
conv_root = .true.
call cgsolve_all (ch_psi_all, cg_psi, et(1,ikmk), dvpsi, dpsi, &
h_diag, npwx, npwq, thresh, ik, lter, conv_root, &
anorm, nbnd_occ(ikk), npol )
IF (isolv==2) THEN
vrs(:,2:4)=-vrs(:,2:4)
IF (okvan) deeq_nc(:,:,:,:)=deeq_nc_save(:,:,:,:,1)
ENDIF
ltaver = ltaver + lter
lintercall = lintercall + 1
if (.not.conv_root) WRITE( stdout, '(5x,"kpoint",i4," ibnd",i4, &
& " solve_linter: root not converged ",es10.3)') &
& ik , ibnd, anorm
!
! writes delta_psi on iunit iudwf, k=kpoint,
!
! if (nksq.gt.1 .or. npert(irr).gt.1)
call save_buffer (dpsi, lrdwf, iudwf, nrec)
!
! calculates dvscf, sum over k => dvscf_q_ipert
!
weight = wk (ikk)
IF (nsolv==2) weight=weight/2.0_DP
IF (noncolin) THEN
call incdrhoscf_nc(drhoscf(1,1,ipert),weight,ik, &
dbecsum_nc(1,1,1,1,ipert,isolv), dpsi, rsign)
ELSE
call incdrhoscf (drhoscf(1,current_spin,ipert), weight, ik, &
dbecsum(1,1,current_spin,ipert), dpsi)
END IF
! on perturbations
enddo
! on isolv
END DO
! on k-points
enddo
vrs(:, 2:4) = -vrs(:, 2:4)
IF (okvan) deeq_nc(:,:,:,:) = deeq_nc_save(:,:,:,:,2)
ENDIF
!
! set threshold for iterative solution of the linear system
!
IF (first_iter) THEN
thresh = 1.0d-2
ELSE
thresh = min (1.d-1 * sqrt (dr2), 1.d-2)
ENDIF
!
! Compute drhoscf, the charge density response to the total potential
!
CALL sternheimer_kernel(first_iter, isolv==2, npe, lrbar, iubar, &
thresh, dvscfins, all_conv, averlt, drhoscf, dbecsum, &
dbecsum_nc(:,:,:,:,:,isolv))
!
! reset the original magnetic field if it was changed
!
IF (isolv == 2) THEN
IF (.NOT. first_iter) THEN
dvscfins(:, 2:4, :) = -dvscfins(:, 2:4, :)
IF (okvan) int3_nc(:,:,:,:,:) = int3_save(:,:,:,:,:,1)
ENDIF
vrs(:, 2:4) = -vrs(:, 2:4)
IF (okvan) deeq_nc(:,:,:,:) = deeq_nc_save(:,:,:,:,1)
ENDIF
!
END DO ! isolv
!
IF (nsolv==2) THEN
drhoscf = drhoscf / 2.0_DP
dbecsum = dbecsum / 2.0_DP
dbecsum_nc = dbecsum_nc / 2.0_DP
ENDIF
!
! The calculation of dbecsum is distributed across processors (see addusdbec)
! Sum over processors the contributions coming from each slice of bands
@ -635,14 +530,7 @@ SUBROUTINE solve_linter (irr, imode0, npe, drhoscf)
int3_nc(:,:,:,:,:)=int3_save(:,:,:,:,:,1)
ENDIF
END IF
#if defined(__MPI)
aux_avg (1) = DBLE (ltaver)
aux_avg (2) = DBLE (lintercall)
call mp_sum ( aux_avg, inter_pool_comm )
averlt = aux_avg (1) / aux_avg (2)
#else
averlt = DBLE (ltaver) / lintercall
#endif
!
tcpu = get_clock ('PHONON')
WRITE( stdout, '(/,5x," iter # ",i3," total cpu time :",f8.1, &
@ -694,7 +582,6 @@ SUBROUTINE solve_linter (irr, imode0, npe, drhoscf)
CALL apply_dpot_deallocate()
if (allocated(ldoss)) deallocate (ldoss)
if (allocated(ldos)) deallocate (ldos)
deallocate (h_diag)
deallocate (dbecsum)
IF (okpaw) THEN
if (allocated(becsum1)) deallocate (becsum1)

3
PHonon/PH/zstar_eu.f90
View File

@ -28,8 +28,7 @@ subroutine zstar_eu
USE qpoint, ONLY : npwq, nksq, ikks
USE eqv, ONLY : dvpsi, dpsi
USE efield_mod, ONLY : zstareu0, zstareu
USE units_ph, ONLY : iudwf, lrdwf
USE units_lr, ONLY : iuwfc, lrwfc
USE units_lr, ONLY : iuwfc, lrwfc, iudwf, lrdwf
USE control_lr, ONLY : nbnd_occ
USE control_ph, ONLY : done_zeu
USE ph_restart, ONLY : ph_writefile

6
PHonon/PH/zstar_eu_us.f90
View File

@ -33,9 +33,9 @@ subroutine zstar_eu_us
USE efield_mod, ONLY : zstareu0
USE phus, ONLY : becsumort
USE modes, ONLY : u, npert, nirr
USE units_ph, ONLY : lrdwf, iucom, lrcom, lrebar, iuebar, &
lrdrhous, iudrhous, iudwf
USE units_lr, ONLY : iuwfc, lrwfc
USE units_ph, ONLY : iucom, lrcom, lrebar, iuebar, &
lrdrhous, iudrhous
USE units_lr, ONLY : iuwfc, lrwfc, iudwf, lrdwf
USE mp_pools, ONLY : nproc_pool, npool
USE control_lr, ONLY : nbnd_occ
USE lrus, ONLY : int3, int3_paw

118
PP/src/Makefile
View File

@ -52,10 +52,8 @@ write_io_header.o \
write_hamiltonians.o \
xc_vdW_scale_mod.o # added by Yang Jiao
PWOBJS = ../../PW/src/libpw.a
QEMODS = ../../Modules/libqemod.a ../../upflib/libupf.a ../../KS_Solvers/libks_solvers.a \
../../FFTXlib/libqefft.a ../../LAXlib/libqela.a ../../UtilXlib/libutil.a \
../../dft-d3/libdftd3qe.a ../../XClib/xc_lib.a
PWOBJS = ../../PW/src/libpw.a ../../KS_Solvers/libks_solvers.a ../../dft-d3/libdftd3qe.a
QEMODS = $(BASEMODS)
MODULES = $(PWOBJS) $(QEMODS)
@ -73,153 +71,153 @@ libpp.a : $(PPOBJS)
$(RANLIB) $@
qexml.x: qexml.o qexml_example.o
$(LD) $(LDFLAGS) -o $@ qexml_example.o qexml.o $(QEMODS) $(LIBOBJS)
$(LD) $(LDFLAGS) -o $@ qexml_example.o qexml.o $(QEMODS)
open_grid.x : open_grid.o libpp.a $(MODULES) $(LIBOBJS)
open_grid.x : open_grid.o libpp.a $(MODULES)
$(LD) $(LDFLAGS) -o $@ \
open_grid.o libpp.a $(MODULES) $(LIBOBJS) $(QELIBS)
open_grid.o libpp.a $(MODULES) $(QELIBS)
- ( cd ../../bin ; ln -fs ../PP/src/$@ . )
average.x : average.o libpp.a $(MODULES) $(LIBOBJS)
average.x : average.o libpp.a $(MODULES)
$(LD) $(LDFLAGS) -o $@ \
average.o libpp.a $(MODULES) $(LIBOBJS) $(QELIBS)
average.o libpp.a $(MODULES) $(QELIBS)
- ( cd ../../bin ; ln -fs ../PP/src/$@ . )
bands.x : bands.o libpp.a $(MODULES) $(LIBOBJS)
bands.x : bands.o libpp.a $(MODULES)
$(LD) $(LDFLAGS) -o $@ \
bands.o libpp.a $(MODULES) $(LIBOBJS) $(QELIBS)
bands.o libpp.a $(MODULES) $(QELIBS)
- ( cd ../../bin ; ln -fs ../PP/src/$@ . )
dos.x : dos.o libpp.a $(MODULES) $(LIBOBJS)
dos.x : dos.o libpp.a $(MODULES)
$(LD) $(LDFLAGS) -o $@ \
dos.o libpp.a $(MODULES) $(LIBOBJS) $(QELIBS)
dos.o libpp.a $(MODULES) $(QELIBS)
- ( cd ../../bin ; ln -fs ../PP/src/$@ . )
pawplot.x : pawplot.o libpp.a libpp.a $(MODULES) $(LIBOBJS)
pawplot.x : pawplot.o libpp.a libpp.a $(MODULES)
$(LD) $(LDFLAGS) -o $@ \
pawplot.o libpp.a $(MODULES) $(LIBOBJS) $(QELIBS)
pawplot.o libpp.a $(MODULES) $(QELIBS)
plan_avg.x : plan_avg.o libpp.a $(MODULES) $(LIBOBJS)
plan_avg.x : plan_avg.o libpp.a $(MODULES)
$(LD) $(LDFLAGS) -o $@ \
plan_avg.o libpp.a $(MODULES) \
$(LIBOBJS) $(QELIBS)
$(QELIBS)
- ( cd ../../bin ; ln -fs ../PP/src/$@ . )
plotband.x : plotband.o libpp.a $(MODULES) $(LIBOBJS)
plotband.x : plotband.o libpp.a $(MODULES)
$(LD) $(LDFLAGS) -o $@ \
plotband.o libpp.a $(MODULES) $(LIBOBJS) $(QELIBS)
plotband.o libpp.a $(MODULES) $(QELIBS)
- ( cd ../../bin ; ln -fs ../PP/src/$@ . )
plotproj.x : plotproj.o libpp.a $(MODULES) $(LIBOBJS)
plotproj.x : plotproj.o libpp.a $(MODULES)
$(LD) $(LDFLAGS) -o $@ \
plotproj.o libpp.a $(MODULES) $(LIBOBJS) $(QELIBS)
plotproj.o libpp.a $(MODULES) $(QELIBS)
- ( cd ../../bin ; ln -fs ../PP/src/$@ . )
plotrho.x : plotrho.o libpp.a $(MODULES) $(LIBOBJS)
plotrho.x : plotrho.o libpp.a $(MODULES)
$(LD) $(LDFLAGS) -o $@ \
plotrho.o libpp.a $(MODULES) $(LIBOBJS) $(QELIBS)
plotrho.o libpp.a $(MODULES) $(QELIBS)
- ( cd ../../bin ; ln -fs ../PP/src/$@ . )
pmw.x : poormanwannier.o libpp.a $(MODULES) $(LIBOBJS)
pmw.x : poormanwannier.o libpp.a $(MODULES)
$(LD) $(LDFLAGS) -o $@ \
poormanwannier.o libpp.a $(MODULES) \
$(LIBOBJS) $(QELIBS)
$(QELIBS)
- ( cd ../../bin ; ln -fs ../PP/src/$@ . )
pp.x : postproc.o libpp.a $(MODULES) $(LIBOBJS)
pp.x : postproc.o libpp.a $(MODULES)
$(LD) $(LDFLAGS) -o $@ \
postproc.o libpp.a $(MODULES) \
$(LIBOBJS) $(QELIBS)
$(QELIBS)
- ( cd ../../bin ; ln -fs ../PP/src/$@ . )
projwfc.x : projwfc.o libpp.a $(MODULES) $(LIBOBJS)
projwfc.x : projwfc.o libpp.a $(MODULES)
$(LD) $(LDFLAGS) -o $@ \
projwfc.o libpp.a $(MODULES) $(LIBOBJS) $(QELIBS)
projwfc.o libpp.a $(MODULES) $(QELIBS)
- ( cd ../../bin ; ln -fs ../PP/src/$@ . )
pw2wannier90.x : pw2wannier90.o libpp.a $(MODULES) $(LIBOBJS)
pw2wannier90.x : pw2wannier90.o libpp.a $(MODULES)
$(LD) $(LDFLAGS) -o $@ \
pw2wannier90.o libpp.a $(MODULES) $(LIBOBJS) $(QELIBS)
pw2wannier90.o libpp.a $(MODULES) $(QELIBS)
- ( cd ../../bin ; ln -fs ../PP/src/$@ . )
pw2critic.x : pw2critic.o libpp.a $(MODULES) $(LIBOBJS)
pw2critic.x : pw2critic.o libpp.a $(MODULES)
$(LD) $(LDFLAGS) -o $@ \
pw2critic.o libpp.a $(MODULES) $(LIBOBJS) $(QELIBS)
pw2critic.o libpp.a $(MODULES) $(QELIBS)
- ( cd ../../bin ; ln -fs ../PP/src/$@ . )
wfck2r.x : wfck2r.o libpp.a $(MODULES) $(LIBOBJS)
wfck2r.x : wfck2r.o libpp.a $(MODULES)
$(LD) $(LDFLAGS) -o $@ \
wfck2r.o libpp.a $(MODULES) $(LIBOBJS) $(QELIBS)
wfck2r.o libpp.a $(MODULES) $(QELIBS)
- ( cd ../../bin ; ln -fs ../PP/src/$@ . )
initial_state.x : initial_state.o libpp.a $(MODULES) $(LIBOBJS)
initial_state.x : initial_state.o libpp.a $(MODULES)
$(LD) $(LDFLAGS) -o $@ \
initial_state.o libpp.a $(MODULES) \
$(LIBOBJS) $(QELIBS)
$(QELIBS)
- ( cd ../../bin ; ln -fs ../PP/src/$@ . )
pw2gw.x : pw2gw.o libpp.a $(MODULES) $(LIBOBJS)
pw2gw.x : pw2gw.o libpp.a $(MODULES)
$(LD) $(LDFLAGS) -o $@ \
pw2gw.o libpp.a $(MODULES) $(LIBOBJS) $(QELIBS)
pw2gw.o libpp.a $(MODULES) $(QELIBS)
- ( cd ../../bin ; ln -fs ../PP/src/$@ . )
sumpdos.x : sumpdos.o
$(LD) $(LDFLAGS) -o $@ sumpdos.o $(MODULES) $(LIBOBJS) $(QELIBS)
$(LD) $(LDFLAGS) -o $@ sumpdos.o $(MODULES) $(QELIBS)
- ( cd ../../bin ; ln -fs ../PP/src/$@ . )
epsilon.x : epsilon.o libpp.a $(MODULES) $(LIBOBJS)
epsilon.x : epsilon.o libpp.a $(MODULES)
$(LD) $(LDFLAGS) -o $@ epsilon.o libpp.a $(MODULES) \
$(LIBOBJS) $(QELIBS)
$(QELIBS)
- ( cd ../../bin ; ln -fs ../PP/src/$@ . )
wannier_ham.x : wannier_ham.o libpp.a $(MODULES) $(LIBOBJS)
wannier_ham.x : wannier_ham.o libpp.a $(MODULES)
$(LD) $(LDFLAGS) -o $@ \
wannier_ham.o libpp.a $(MODULES) $(LIBOBJS) $(QELIBS)
wannier_ham.o libpp.a $(MODULES) $(QELIBS)
- ( cd ../../bin ; ln -fs ../PP/src/$@ . )
wannier_plot.x : wannier_plot.o libpp.a $(MODULES) $(LIBOBJS)
wannier_plot.x : wannier_plot.o libpp.a $(MODULES)
$(LD) $(LDFLAGS) -o $@ \
wannier_plot.o libpp.a $(MODULES) $(LIBOBJS) $(QELIBS)
wannier_plot.o libpp.a $(MODULES) $(QELIBS)
- ( cd ../../bin ; ln -fs ../PP/src/$@ . )
molecularpdos.x : molecularpdos.o libpp.a $(MODULES) $(LIBOBJS)
molecularpdos.x : molecularpdos.o libpp.a $(MODULES)
$(LD) $(LDFLAGS) -o $@ \
molecularpdos.o libpp.a $(MODULES) $(LIBOBJS) $(QELIBS)
molecularpdos.o libpp.a $(MODULES) $(QELIBS)
- ( cd ../../bin ; ln -fs ../PP/src/$@ . )
pw2bgw.x : pw2bgw.o libpp.a $(MODULES) $(LIBOBJS)
pw2bgw.x : pw2bgw.o libpp.a $(MODULES)
$(LD) $(LDFLAGS) -o $@ \
pw2bgw.o libpp.a $(MODULES) $(LIBOBJS) $(QELIBS)
pw2bgw.o libpp.a $(MODULES) $(QELIBS)
- ( cd ../../bin ; ln -fs ../PP/src/$@ . )
#bgw2pw.x : bgw2pw.o libpp.a $(MODULES) $(LIBOBJS)
#bgw2pw.x : bgw2pw.o libpp.a $(MODULES)
# $(LD) $(LDFLAGS) -o $@ \
# bgw2pw.o libpp.a $(MODULES) $(LIBOBJS) $(QELIBS)
# bgw2pw.o libpp.a $(MODULES) $(QELIBS)
# - ( cd ../../bin ; ln -fs ../PP/src/$@ . )
fs.x : fermisurface.o libpp.a $(MODULES) $(LIBOBJS)
fs.x : fermisurface.o libpp.a $(MODULES)
$(LD) $(LDFLAGS) -o $@ \
fermisurface.o libpp.a $(MODULES) $(LIBOBJS) $(QELIBS)
fermisurface.o libpp.a $(MODULES) $(QELIBS)
- ( cd ../../bin ; ln -fs ../PP/src/$@ . )
tldeps :
if test -n "$(TLDEPS)" ; then \
( cd ../.. ; $(MAKE) $(TLDEPS) || exit 1 ) ; fi
fermi_velocity.x : fermi_velocity.o libpp.a $(MODULES) $(LIBOBJS)
fermi_velocity.x : fermi_velocity.o libpp.a $(MODULES)
$(LD) $(LDFLAGS) -o $@ \
fermi_velocity.o libpp.a $(MODULES) $(LIBOBJS) $(QELIBS)
fermi_velocity.o libpp.a $(MODULES) $(QELIBS)
- ( cd ../../bin ; ln -fs ../PP/src/$@ . )
fermi_proj.x : fermi_proj.o libpp.a $(MODULES) $(LIBOBJS)
fermi_proj.x : fermi_proj.o libpp.a $(MODULES)
$(LD) $(LDFLAGS) -o $@ \
fermi_proj.o libpp.a $(MODULES) $(LIBOBJS) $(QELIBS)
fermi_proj.o libpp.a $(MODULES) $(QELIBS)
- ( cd ../../bin ; ln -fs ../PP/src/$@ . )
# Added by Yang Jiao
ppacf.x : ppacf.o libpp.a $(MODULES) $(LIBOBJS)
ppacf.x : ppacf.o libpp.a $(MODULES)
$(LD) $(LDFLAGS) -o $@ \
ppacf.o libpp.a $(MODULES) $(LIBOBJS) $(QELIBS)
ppacf.o libpp.a $(MODULES) $(QELIBS)
- ( cd ../../bin ; ln -fs ../PP/src/$@ . )
clean :

Some files were not shown because too many files have changed in this diff Show More