Piezoelectric tensor calculation fail

[gutubu]

Dear all:

I am trying to do a piezoelectric tensor calculation following the tutorial elastic - abinit using Abinit 10.4.5.. Now in step 2 (calculation of several 2nd order derivatives of the total energy) it seems Ok until the perturbation reaches the last atom and the calculation stopped without outputting any DDB file? Please see the pasted .abi and .log files. Similar situation occured twice after a one-week run of the program. Thank you.

#KGN step 2

ndtset 3

###############################################################################

Set 1 : Initial self-consistent run

###############################################################################
kptopt1 1
tolvrs1 1.0d-18 #need excellent convergence of GS quantities for RF runs
########################################################################

Set 2 : Calculate the ddk wf’s - needed for piezoelectric tensor and Born effective charges in dataset 3

#######################################################################
getwfk2 -1
iscf2 -3 #this option is needed for ddk
kptopt2 2 #use time-reversal symmetry only for k points
nqpt2 1 #one wave vector will be specified
qpt2 0 0 0 #need to specify gamma point
rfelfd2 2 #set for ddk wf’s only
tolwfr2 1.0d-20 #only wf convergence can be monitored here
#############################################################################

Set 3 : response-function calculations for all needed perturbations

############################################################################
getwfk3 -2
getddk3 -1
kptopt3 2 #use time-reversal symmetry only for k points
nqpt3 1
qpt3 0 0 0
rfphon3 1 #do atomic displacement perturbation
rfstrs3 3 #do strain perturbation
tolvrs3 1.0d-10 #need reasonable convergence of 1st-order quantities
##################################################################

Common input data

##################################################################
acell 1.9630015360E+01 1.8731783186E+01 3.0242540769E+01 Bohr

rprim 9.9998997866E-01 0.0000000000E+00 4.4768932538E-03
0.0000000000E+00 1.0000000000E+00 0.0000000000E+00
-1.7400338913E-01 0.0000000000E+00 9.8474505359E-01

#Definition of the atom types and atoms
ntypat 4
znucl 8.00000 7.00000 6.00000 1.00000
natom 156
typat 1 1 1 1 1 1 1 1 1 1 1 1 2 2 2 2 3 3 3 3
3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3
3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3
3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3
3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 4 4 4 4
4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4
4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4
4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4

xred ####################COPY RELAXED RESULT FROM PREVIOUS CALCULATION
2.0242339844E-01 3.0532617535E-01 5.0137591005E-01
7.9757660156E-01 8.0532617535E-01 4.9862408995E-01

……..

#Gives the number of bands, explicitely (do not take the default)
nband 248 # For an insulator (if described correctly as an

insulator by DFT), conduction bands should not

be included in response-function calculations

#Definition of the plane wave basis set
ecut 22.0 # Maximum kinetic energy cutoff (Hartree)
ecutsm 0.5 # Smoothing energy needed for lattice paramete

optimization. This will be retained for

consistency throughout.

#Definition of the k-point grid
kptopt 1 # Use symmetry and treat only inequivalent points
ngkpt 1 1 1 # 1x1x1 Monkhorst-Pack grid (gamma)
nshiftk 1 # Use one copy of grid only (default)
shiftk 0.0 0.0 0.0 # This choice of origin for the k point grid

preserves the hexagonal symmetry of the grid,

which would be broken by the default choice.

#Definition of the self-consistency procedure
diemac 3.0 # Model dielectric preconditioner (model DIElectric MACroscopic constant)
#, For wider gap insulators, use 2.0 … 4.0
nstep 600 # Maxiumum number of SCF iterations

enforce calculation of forces at each SCF step

optforces 1

pp_dirpath “$ABI_PSPDIR”
pseudos “O.psp8, N.psp8, C.psp8, H.psp8”

KGN_elast2_1_248band.log (2.8 MB)

[mverstra]

Hi,

You get a segmentation fault, could be a memory issue or configuration of your mpi/ compiler/ cluster: it complains the allocation of a chunk is not aligned with the rest of the memory. This happens during the non stationary calculation of off diagonal d2E/dRdE derivatives. Overall your calculation looks small so it’s probably not a lack of memory.

The most efficient thing to do is run one perturbation at a time, as separate jobs, and then merge at the end. This avoids loading all perturbative in memory at once, and writes 1ddb file at a time so nothing is lost. If you use abipy it does this for you automatically. The main issue is identifying the symmetry irreducible perturbations. You can look through your log file. To start you can run just the last perturbation, which failed, to see what happens.

It’s probable that if you run all pert separately they will complete, and you will be able to proceed. NB: you only have 1 kpt, this is probably unconverged

[gutubu]

Dear professor：

Thank you for your prompt reply. I am very new to abinit and am not using abipy at the stage.
Would you kindly see the attached .abi file and instruct on how to set properly the variable(s) for “run one perturbation at a time”? Also, my previous run generated quite a lot files, including 1WFs, DENs, POTs, is there a way to make use of these files so that a re-calculation can start reading these files and thus accelerated?

ZX
KGN_elast2_1_adjusted.abi (13.0 KB)

[mverstra]

1. parse the log file to get the list of perturbations which are irreducible by symmetry

2. make 1 directory per perturbation, copy the input file, and make symbolic links to the psp files and WFK ground state files (+ the DDK wave functions if you are doing DDE perturbations)

3. in each subdirectory edit input file to only run the DFPT calculation and set

ndtset 1

jdtset 3

getwfk3 1

rfatpol3 i i # here i is the atom being perturbed according to the list in 1)

rfdir3 0 1 0 # example for direction 2 = reduced y axis. Adjust according to list in 1)

rfelfd3 3 # in case of DDE perturbation at Gamma

qpt3 qx qy qz # if you are running at finite q, set it here - probably not your case for just Gamma

4. run each of the calculations separately (try 1 first as a test before launching everything else)

5. collect all of the DDB files and merge

[gutubu]

Dear professor:

Thank you for your reply.
I am trying to use abipy for the automatic flow and now I am in Dataset 1. I found that

================================================================================
== DATASET 1 ==================================================================

• mpi_nproc: 18, omp_nthreads: -1 (-1 if OMP is not activated)

I’m only using 18 cores for this calculation, which I have no clue where I have set that. Perhaps from

acell 1.0000000000E+00 1.0000000000E+00 1.0000000000E+00 Bohr
amu 1.59994000E+01 1.40067400E+01 1.20110000E+01
1.00794000E+00
diemac 3.00000000E+00
ecut 2.20000000E+01 Hartree
ecutsm 5.00000000E-01 Hartree
fftalg 312
istwfk 2 0 3 0 0 0 6 0 7 4
0 5 0 0 0 8 0 9
ixc 11
kpt 0.00000000E+00 0.00000000E+00 0.00000000E+00
2.50000000E-01 0.00000000E+00 0.00000000E+00
5.00000000E-01 0.00000000E+00 0.00000000E+00
0.00000000E+00 2.50000000E-01 0.00000000E+00
2.50000000E-01 2.50000000E-01 0.00000000E+00
5.00000000E-01 2.50000000E-01 0.00000000E+00
0.00000000E+00 5.00000000E-01 0.00000000E+00
2.50000000E-01 5.00000000E-01 0.00000000E+00
5.00000000E-01 5.00000000E-01 0.00000000E+00
0.00000000E+00 0.00000000E+00 5.00000000E-01
2.50000000E-01 0.00000000E+00 5.00000000E-01
5.00000000E-01 0.00000000E+00 5.00000000E-01
0.00000000E+00 2.50000000E-01 5.00000000E-01
2.50000000E-01 2.50000000E-01 5.00000000E-01
5.00000000E-01 2.50000000E-01 5.00000000E-01
0.00000000E+00 5.00000000E-01 5.00000000E-01
2.50000000E-01 5.00000000E-01 5.00000000E-01
5.00000000E-01 5.00000000E-01 5.00000000E-01
kptrlatt 4 0 0 0 4 0 0 0 2
kptrlen 6.04850815E+01
P mkmem 1
natom 156
nband 248
ngfft 90 80 128
nkpt 18
nstep 600
nsym 2
ntypat 4
…
the program automatically set the nproc it would use? Is there any way to modify this or is it not necessary? Thank you.

ZX

[mverstra]

I guess it will try to optimize with not. You can force it with *npkpt* in the input file, or/and by running it with

mpirun –nproc 9

Or whatever number you like, preferably a divisior of nkpt.

If you use more than nkpt it will add band parallelization, be sure to give it 18*N processors, where N is a divisior of nband. To make this easier, you can increase nband in all input files to something round, adding some unoccupied bands.

[gutubu]

Dear professor:

I ended up setting “max_ncores_used” to 79 which is a divisor of my nband of 237, and the flow is running smoothly so far with 79 cores with the “max_njobs_inqueue” set to 1 in the scheduler.yml..Thank you for your previous answers.

ZX

[gutubu]

Dear professor:

How should I set the tolwfr in the DDK calculation? In the run_elastic.py I set it like this

which seemed not read and passed to the .in file for any task(neither scf nor ddk).

Here is the .py script. Would you kindly help with this? Thank you.

#!/usr/bin/env python
r"“”
Elastic constants and piezoelectric tensor with DFPT

This example shows how to use AbiPy to calculate physical properties
related to strain for an insulator.

- the rigid-atom elastic tensor
- the rigid-atom piezoelectric tensor (insulators only)
- the internal strain tensor
- the atomic relaxation corrections to the elastic and piezoelectric tensor

Here we follow the discussion presented in
in the the official tutorial <https://docs.abinit.org/tutorial/elastic/>_

The DDB file with all the perturbations will be produced automatically at the end of the run
and saved in flow_elastic/w0/outdata/out_DDB.
“”"
import sys
import os
import numpy as np
import abipy.abilab as abilab
import abipy.data as abidata
from abipy import flowtk
from abipy.flowtk.launcher import PyLauncher

def make_scf_input():
“”"
This function constructs the input file for the GS calculation of
KGN in P21 spg.
In principle, the stucture should be relaxed before starting the calculation
“”"

# Initialize structure. Use enough significant digits
# so that Abinit will recognize the correct spacegroup
# (Hexagonal and rhombohedral lattices are a bit problematic).
structure = abilab.Structure.from_abivars(
    acell=[1.9630015360E+01, 1.8731783186E+01, 3.0242540769E+01],
    natom=156,
    ntypat=4,
    rprim=[
        9.9998997866E-01, 0.0000000000E+00, 4.4768932538E-03,
        0.0000000000E+00, 1.0000000000E+00, 0.0000000000E+00,
        -1.7400338913E-01, 0.0000000000E+00, 9.8474505359E-01
    ],
    typat=[
        1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 2, 2, 2, 2, 3, 3, 3, 3,
        3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3,
        3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3,
        3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3,
        3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 4, 4, 4, 4,
        4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4,
        4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4,
        4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4
    ],
    xred=[
        2.0242339844E-01, 3.0532617535E-01, 5.0137591005E-01,
        7.9757660156E-01, 8.0532617535E-01, 4.9862408995E-01,
        3.0218345846E-01, 8.4767988126E-01, 4.9554416146E-01,
        ...
    ],
    znucl=[8, 7, 6, 1],
)

#pseudos = abidata.pseudos("8o.pspnc", "1h.paw", "7n.pspnc", "6c.pspnc")
pseudos = abidata.pseudos("O.psp8", "H.psp8", "N.psp8", "C.psp8")
gs_inp = abilab.AbinitInput(structure, pseudos=pseudos)

# Set other important variables (consistent with tutorial)
# All the other DFPT runs will inherit these parameters.
gs_inp.set_vars(
    nband=237, # Number of bands (should be checked carefully)
    ecut=20.22, # Energy cutoff
    ecutsm=0.5,        # Important when performing structural optimization
    # with variable cell. All DFPT calculations should use
    # the same value to be consistent.
    ngkpt=[1, 1, 1], # k-mesh (to be converged)
    nshiftk=1, # No k-shift, Gamma-centered grid
    shiftk=[0.0, 0.0, 0.0],   # Gamma-centered k-mesh
    diemac=3.0, # Dielectric constant for Coulomb cutoff
    nstep=1200, # Maximum number of SCF cycles
    tolvrs=1.0e-18, # Very tight tolerance on residuals
    autoparal=1,      # Automatically determine the best way to parallelize
)
return gs_inp

def build_flow(options):
“”"
Create a Flow for phonon calculations. The flow has one work with:

    - 1 GS Task
    - 3 DDK Task
    - 4 Phonon Tasks (Gamma point)
    - 6 Elastic tasks (3 uniaxial + 3 shear strain)

The Phonon tasks and the elastic task will read the DDK produced at the beginning
"""
# Working directory (default is the name of the script with '.py' removed and "run_" replaced by "flow_")
if not options.workdir:
    options.workdir = os.path.basename(sys.argv[0]).replace(".py", "").replace("run_", "flow_")

flow = flowtk.Flow(workdir=options.workdir)

# Build input for GS calculation and register the first work.
scf_input = make_scf_input()

elast_work = flowtk.ElasticWork.from_scf_input(
    scf_input, 
    with_relaxed_ion=True, 
    with_piezo=True, 
    tolerances={"tolwfr": 1.0e-20}   # NSCF tolerance
)
flow.register_work(elast_work)

return flow

This block generates the thumbnails in the AbiPy gallery.

You can safely REMOVE this part if you are using this script for production runs.

if os.getenv(“READTHEDOCS”, False):
name = None
import tempfile
options = flowtk.build_flow_main_parser().parse_args([“-w”, tempfile.mkdtemp()])
build_flow(options).graphviz_imshow()

@flowtk.flow_main
def main(options):
“”"
This is our main function that will be invoked by the script.
flow_main is a decorator implementing the command line interface.
Command line args are stored in options.
“”"
return build_flow(options)

if name == “main”:
sys.exit(main())

[gutubu]

Dear Professor:

My calculation of the piezoelectric tensor using abipy script ended up terminated at the “elasticTask” after about one-week’s run.

The .err of the first “elasticTask” stated exactly the same problem as the last time when I tried the step-by-step calculation:

Please help. I am also uploading the run_elastic.py as a .txt so that if you got a chance please take a look.
run_elastic.txt (15.0 KB)

Thank you.

[mverstra]

Hi again,

you are running a large system with still quite a few k-points, so things are heavy. The run is failing in the non stationary calculation of off-diagonal energy derivatives with E fields, and a lot of stuff is allocated there, so I would bet it’s a memory problem. Apart from recommending you try different parallel distributions or more cores, I don’t know what to suggest else.

If you are already done with the SCF convergence of the DFPT for this elastic part you can in principle restart from there by moving the 1WF files and running with iscf -3 and ird1wf.

Basically, you should try an small/intermediate sized system first, to make sure your workflow is complete and functional.

M.

[gutubu]

Thank you.
I changed to a smaller system with only 38 atoms this time, and used actually fewer cores and the flow completed successfully. However, the result of the piezo tensor is way lower (5 orders of magnitude) than my experimental, which is about 138 pC/N.
====================== piezoelectric tensors available ======================
[PIEZO_RELAXED]
relaxed-ion piezoelectric tensor in Voigt notation (shape: (3, 6))
Units: c/m^2, set to zero below: 1e-05, fit_to_structure: True

      xx        yy        zz        yz        xz        xy

Px -0.000689 -0.012647 -0.019409 0.028695 0.002073 -0.033748
Py -0.023563 -0.070751 0.020530 -0.038228 -0.004476 0.009140
Pz -0.010013 -0.003386 -0.076848 -0.018777 0.010102 -0.010849

[PIEZO_CLAMPED]
clamped-ion piezoelectric tensor in Voigt notation (shape: (3, 6))
Units: c/m^2, set to zero below: 1e-05, fit_to_structure: True

      xx        yy        zz        yz        xz        xy

Px -0.001659 -0.000721 0.001214 0.002138 0.025802 0.027896
Py 0.005149 0.083966 -0.047245 0.011327 0.003351 -0.002782
Pz 0.003288 -0.041556 0.084424 0.000987 -0.001542 0.003442

One thing I noticed is that in the tutorial Elastic constants and piezoelectric tensor with DFPT — abipy 0.9.8 documentation, the “Has (at least one) electric-field perturbation: True” while mine is False. Other than that everything went smoothly.