"Subroutine Unknown" and bad results for elastic calculation

vistawanted · Post by **vistawanted** » Fri Aug 03, 2012 4:34 am

Dear all:
Hi.
I am trying to calculate elastic constants via response method in Abinit 6.12.3, which compiled by ifort in CentOS 6.3 with MKL (Blas+Lapack+FFTW3), the version of ifort is as:

Code: Select all

ifort -V
Intel(R) Fortran Intel(R) 64 Compiler XE for applications running on Intel(R) 64, Version 12.1.3.293 Build 20120212
Copyright (C) 1985-2012 Intel Corporation.  All rights reserved.
FOR NON-COMMERCIAL USE ONLY

I have tried the example for AlAs, which gives reasonable results. But, when I calculated TiC (Fm-3m, or, NaCl structure. Conductor), the warning occurred as in the DATASET 2, corresponding to the DDK calculation:

Code: Select all

 ETOT  5  -698.43961312179     -4.801E-04 1.939E-02 0.000E+00

 Subroutine Unknown:0:WARNING
   New trial energy at line       4 =  -1.142031E+03
     is higher than former: -1.142031E+03

Such a warning kept appearing, although the calculation can be finished. But, after the calculation, I cannot get the elastic properties from anaddb analyzing, the constants are shown as:

Code: Select all

.Version 6.12.3 of ANADDB 
.(MPI version, prepared for a x86_64_linux_intel12.1 computer) 
......
 Elastic Tensor (clamped ion) (unit:10^2GP):

************************************   0.0022444  -0.0003354  -0.0059745
************************************   0.0058939  -0.0000122  -0.0059745
************************************   0.0058940  -0.0003354  -0.0021273
   0.0000000   0.0000000   0.0000000************   0.0000000   0.0000000
   0.0000000   0.0000000   0.0000000   0.0000000************   0.0000000
   0.0000000   0.0000000   0.0000000   0.0000000   0.0000000************

My input file is:

Code: Select all

ndtset  3

# GS calculation
iscf1  7
kptopt1  1
tolvrs1  1.0d-18

# DDK Calculation for Wavefunctions
getwfk2  -1
iscf2  -3
kptopt2  2
nqpt2  1
qpt2  0  0  0
rfelfd2  2
rfdir2  1  1  1
tolwfr2  1.0d-20

# Linear Lesponse for Perturabtions
getddk3  -1
getwfk3  -2
iscf3  7
kptopt3  2
nqpt3  1
qpt3  0  0  0
rfphon3  1
rfatpol3  1  2
rfstrs3  3
rfdir3  1  1  1
tolvrs3  1.0d-10

# General Parameters
acell  3*8.254628275
ntypat  2
znucl  22  6
natom  2
typat  1  2
rprim
   0.0  0.5  0.5
   0.5  0.0  0.5
   0.5  0.5  0.0
xred
   0.0  0.0  0.0
   0.5  0.5  0.5
nband  13
ecut  45
ecutsm  0.01
occopt  4
tsmear  0.01
kptopt  1
ngkpt  13  13  13
nshiftk  4
shiftk
   0.5  0.0  0.0
   0.0  0.5  0.0
   0.0  0.0  0.5
   0.5  0.5  0.5
nstep  100
iscf  7
diemac  5.0
optforces  1

The pseudo-potentials I was used is, Ti:

Code: Select all

Goedecker pseudopotential for Ti
   22  12  070301 zatom,zion,pspdat
10 11  2 2 2001 0  pspcod,pspxc,lmax,lloc,mmax,r2well
     0.38000000    2     8.71144218    -0.70028677                                  rloc nloc c1 c2
    3                                                                               nnonloc
     0.33777078    2     2.57526386     3.69297065                                  rs ns hs11 hs12
                                       -4.76760461                                             hs22
     0.24253135    2    -4.63054123     8.87087502                                  rp np hp11 hp12
                                      -10.49616087                                             hp22
                         0.20121752    -0.06088315                                        kp11 kp12
                                        0.07203791                                             kp22
     0.24331694    1    -9.40665268                                                 rd nd hd11
                         0.00584456                                                       kd11

carbon:

Code: Select all

Goedecker pseudopotential for C
    6   4  070301 zatom,zion,pspdat
10 11  1 2 2001 0  pspcod,pspxc,lmax,lloc,mmax,r2well
     0.33847124    2    -8.80367398     1.33921085                                  rloc nloc c1 c2
    2                                                                               nnonloc
     0.30257575    1     9.62248665                                                 rs ns hs11
     0.29150694    1     0.00000000                                                 rp np hp11
                         0.00207319                                                       kp11

Would any one kindly give me some hint to solve this problem?
Thanks a lot.

jzwanzig · Post by **jzwanzig** » Fri Aug 03, 2012 3:40 pm

Take out the DDK perturbation (rfelfd 2) and the electric field (rfelfd 3) and re-try. Because your system is a conductor, the electric field perturbation is irrelevant and should fail, and the DDK is an auxiliary to this and should also (if I remember correctly) fail for a conductor. Phonons and strain should be ok. One of the tutorials is explicitly for elastic properties of a metal, be sure to read it and follow the recommendations. Also, by the way, I notice you have diemac set to 5.0, this is a preconditioning value that is appropriate for an insulator, not a metal; the default value of diemac is good for a metal so I'd comment out the diemac line (it might help things converge faster).

vistawanted · Post by **vistawanted** » Sat Aug 04, 2012 1:07 pm

Dear jzwanzig:
Thank you for your reply. I am sorry I haven't realized there is an example for metal before your post.
I have tried two ways to find the elastic constants:
1, GS -> Converge unoccupied wfc -> Response Function ( Similar as the example)
2, GS -> Response Function
Both of them provides same results and a bit far from the experimental results:
--------------C11------C12-----C44 (GPa)
Abinit------564-------145------198
Exp---------500-------100------180
I think this is okey since I used low cut-off energy and rough k-mesh.
But, I think it is valuable to point out, these elastic properties data were obtained from the PBE-FHI pseudo-potential. GTH and HGH pseudo-potentials work very bad in this situation - the calculated elastice constants were higher than several thousands ...... :shock:

jzwanzig · Post by **jzwanzig** » Sat Aug 04, 2012 3:30 pm

You really should compare the pseudopotential performance only in *converged* calculations. You say that the PBE calculation is not converged (or only roughly converged)--then the HGH and GTH calculations will be much more poorly converged, because these types of pseudopotentials are much harder than TM pseudopotentials. In other words, a converged ecut might be 35 Ha for TM, and 80 for HGH and GTH. So the latter two pseudopotentials are not "worse" than TM, you are just running them under much less converged conditions.

vistawanted · Post by **vistawanted** » Sun Aug 12, 2012 5:38 am

jzwanzig wrote:You really should compare the pseudopotential performance only in *converged* calculations. You say that the PBE calculation is not converged (or only roughly converged)--then the HGH and GTH calculations will be much more poorly converged, because these types of pseudopotentials are much harder than TM pseudopotentials. In other words, a converged ecut might be 35 Ha for TM, and 80 for HGH and GTH. So the latter two pseudopotentials are not "worse" than TM, you are just running them under much less converged conditions.

Thank you very much, jzwanzig.
I checked the convergence of GTH in PWscf, they were really harder than normal NCPP......
I used other NCPP with reasonable ecut, and finally found the elastic constants with good agreement to experimental results.

ABINIT Discussion Forums

"Subroutine Unknown" and bad results for elastic calculation

"Subroutine Unknown" and bad results for elastic calculation

Re: "Subroutine Unknown" and bad results for elastic calcula

Re: "Subroutine Unknown" and bad results for elastic calcula

Re: "Subroutine Unknown" and bad results for elastic calcula

Re: "Subroutine Unknown" and bad results for elastic calcula