[Wien] in1new problem

Peter Blaha pblaha at theochem.tuwien.ac.at
Mon Jun 6 16:33:37 CEST 2005


Of course I made an error in my script. It did not do what it should, namely 
put a "search" to the deep semicore states (where you put  <==== should always 
be              this! 
 2   -5.911     0.002 CONT 1  <===========

The new attatched script seems to work.


> > I include the modified (untested) script.
> > Let me know, if this suggestion works.
> >  
> It produces the following case.in1 :
> 
> WFFIL        (WFPRI, SUPWF)
>  7.50       10    4 (R-MT*K-MAX; MAX L IN WF, V-NMT
> .74991   7   0      global e-param with N other choices, napw
> 0    0.663     0.000 CONT 1
> 0   -1.254     0.000 CONT 1
> 1    0.720     0.000 CONT 1
> 1   -0.136     0.000 CONT 1
> 2    0.838     0.000 CONT 1
> 2   -5.911     0.000 CONT 1  <===========
> 3    0.885     0.000 CONT 1
> .74991   6   0      global e-param with N other choices, napw
> 0    0.457     0.000 CONT 1
> 0   -0.212     0.000 CONT 1
> 1    0.868     0.000 CONT 1
> 1   -5.806     0.000 CONT 1  <==========
> 2    0.883     0.000 CONT 1
> 2   -1.036     0.000 CONT 1
> K-VECTORS FROM UNIT:4   -9.0       3.0      emin/emax window
> 
> Compare this to the case.in1 found with the old scheme (starting from the
> case.in1 with the original in1new, at the point where it crashed) :
> 
> WFFIL        (WFPRI, SUPWF)
>  7.50       10    4 (R-MT*K-MAX; MAX L IN WF, V-NMT
> .74994   7   0      global e-param with N other choices, napw
> 0    0.663     0.000 CONT 1
> 0   -1.254     0.005 CONT 1
> 1    0.720     0.000 CONT 1
> 1   -0.137     0.005 CONT 1
> 2    0.838     0.000 CONT 1
> 2   -5.910     0.005 CONT 1   <===========
> 3    0.885     0.000 CONT 1
> .74994   6   0      global e-param with N other choices, napw
> 0    0.457     0.000 CONT 1
> 0   -0.213     0.001 CONT 1
> 1    0.868     0.000 CONT 1
> 1   -5.805     0.005 CONT 1   <========
> 2    0.883     0.000 CONT 1
> 2   -1.036     0.005 CONT 1
> K-VECTORS FROM UNIT:4   -9.0       3.0      emin/emax window
> 
> The energies for the deep semi-core states are identical now. Convergence
> looked OK at first sight, however, there is still a QTL-B crash:
> 
> :DIS  :  CHARGE DISTANCE       0.1318234
> :DIS  :  CHARGE DISTANCE       0.1646432
> :DIS  :  CHARGE DISTANCE       0.1545064
> :DIS  :  CHARGE DISTANCE       0.1624324
> :DIS  :  CHARGE DISTANCE       0.1137017
> :DIS  :  CHARGE DISTANCE       0.0394287
> :DIS  :  CHARGE DISTANCE       0.0909304
> :DIS  :  CHARGE DISTANCE       0.0803767
> :DIS  :  CHARGE DISTANCE       0.0639179
> :DIS  :  CHARGE DISTANCE       0.0504621
> :DIS  :  CHARGE DISTANCE       0.0387721
> :DIS  :  CHARGE DISTANCE       0.0353060
> :DIS  :  CHARGE DISTANCE       0.0204696
> :DIS  :  CHARGE DISTANCE       0.0127861
> :DIS  :  CHARGE DISTANCE       0.0037622
> :DIS  :  CHARGE DISTANCE       0.0056102
> 
> Apparantly this new write_in1 does what it should do. But the linearization
> error slowly grows, and still leads to a crash:
> 
>   QTL-B VALUE .EQ.    2.22915   in Band of energy   -5.83529
>   QTL-B VALUE .EQ.    3.68742   in Band of energy   -5.81296
>   QTL-B VALUE .EQ.    4.39352   in Band of energy   -5.79731
>   QTL-B VALUE .EQ.    4.95376   in Band of energy   -5.78438
>   QTL-B VALUE .EQ.    4.74497   in Band of energy   -5.77338
>   QTL-B VALUE .EQ.    3.79768   in Band of energy   -5.76046
>   QTL-B VALUE .EQ.    5.38024   in Band of energy   -5.74967
>   QTL-B VALUE .EQ.    5.85567   in Band of energy   -5.74256
>   QTL-B VALUE .EQ.    5.83326   in Band of energy   -5.73724
>   QTL-B VALUE .EQ.    5.92836   in Band of energy   -5.73198
>   QTL-B VALUE .EQ.    6.18613   in Band of energy   -5.72743
>   QTL-B VALUE .EQ.    6.25986   in Band of energy   -5.72401
>   QTL-B VALUE .EQ.    6.49172   in Band of energy   -5.72114
>   QTL-B VALUE .EQ.    6.65410   in Band of energy   -5.71830
>   QTL-B VALUE .EQ.    6.77625   in Band of energy   -5.71636
>   QTL-B VALUE .EQ.    6.92212   in Band of energy   -5.71440
> (and it exceeds 7.0 in the following iteration)
> 
> It looks like there is still a difference for these deep states...? The
> problematic case.struct is attached, in case you need (Rkmax=7.5, k-mesh=75).
> Symmetry is artificially reduced to tetragonal for later use.
> 
> Title
> B   LATTICE,NONEQUIV.ATOMS:  2 139 I4/mmm
> MODE OF CALC=RELA
>  7.177134  7.177134 10.150000 90.000000 90.000000 90.000000
> ATOM  -1: X=0.00000000 Y=0.00000000 Z=0.00000000
>          MULT= 1          ISPLIT=-2
> La1        NPT=  781  R0=.000010000 RMT= 2.50        Z:57.0
> LOCAL ROT MATRIX:    1.0000000 0.0000000 0.0000000
>                     0.0000000 1.0000000 0.0000000
>                     0.0000000 0.0000000 1.0000000
> ATOM  -2: X=0.00000000 Y=0.00000000 Z=0.50000000
>          MULT= 1          ISPLIT=-2
> Sb1        NPT=  781  R0=.000010000 RMT= 2.40        Z:51.0
> LOCAL ROT MATRIX:    1.0000000 0.0000000 0.0000000
>                     0.0000000 1.0000000 0.0000000
>                     0.0000000 0.0000000 1.0000000
>  16      NUMBER OF SYMMETRY OPERATIONS
> 
> Thanks,
> Stefaan
> 
> _______________________________________________
> Wien mailing list
> Wien at zeus.theochem.tuwien.ac.at
> http://zeus.theochem.tuwien.ac.at/mailman/listinfo/wien
> 


                                      P.Blaha
--------------------------------------------------------------------------
Peter BLAHA, Inst.f. Materials Chemistry, TU Vienna, A-1060 Vienna
Phone: +43-1-58801-15671             FAX: +43-1-58801-15698
Email: blaha at theochem.tuwien.ac.at    WWW: http://info.tuwien.ac.at/theochem/
--------------------------------------------------------------------------
-------------- next part --------------
#!/bin/csh -f
# reads case.scf and case.in1 and writes a new case.in1
# based on an idea of J.Sofo and J.Fuhr
unalias rm
# define charge limit, when an energy parameter should be set
set qlimit=0.05
# define etest, E-l parameters below etest will be "searched" 
# by original approach using de
set etest=-3.

set updn = ""
set cmplx
unset help

while ($#argv)
  switch ($1)
  case -h:
    set help
    shift; breaksw
  case -ql:
    shift; set qlimit=$1
    shift; breaksw
  case -c:
    set cmplx = c
    shift; breaksw
  case -up:
    set updn = up
    shift; breaksw
  case -dn:
    set updn = dn
    shift; breaksw
  default:
    shift; breaksw
  endsw
end
if ($?help) goto help

set file    = `pwd`
set file    = $file:t

if( ! -e $file.in1$cmplx ) then
    echo " the required file    $file.in1$cmplx    does not exist"
    exit(3)
endif
if( ! -e $file.scf2$updn ) then
    echo " the required file    $file.scf2$updn    does not exist"
    exit(3)
endif

head -2 $file.in1$cmplx > $file.in1"$cmplx"new
set lines=2

#determine EF
set ef1 = (`grep :FER $file.scf2$updn | sed -e 's/:FER.*=//'`)
if ($#ef1 == 0 ) then
            echo " FERMI energy not found in   $file.scf2$updn"
            exit(3)
endif
set eferm=`echo "scale=4;$ef1 - 0.3 " |bc -l`


#loop over atoms 
set natom=`head -2 $file.struct | tail -1 | cut -c28-30`
set iatom=1
while ($iatom <= $natom)
@ lines ++
set nn=`head -$lines $file.in1$cmplx |tail -1`
#find global APW parameter
set apw0=0
if ( $nn[3] == 1 ) set apw0=1
set apws=1
set apwp=1
set apwd=1
set apwf=1
set i0=1

#find l-specific APW parameters from case.in1
while ($i0 <= $nn[2] )
@ lines ++
  set nnn=`head -$lines $file.in1$cmplx |tail -1`
if ($#nnn >= 5 ) then
   if( $nnn[1] == 0 ) set apws=$nnn[5]
   if( $nnn[1] == 1 ) set apwp=$nnn[5]
   if( $nnn[1] == 2 ) set apwd=$nnn[5]
   if( $nnn[1] == 3 ) set apwf=$nnn[5]
endif
@ i0 ++
end

#find energy parameter from case.scf2  (EPL/EPH)
  if ($iatom <= 9) then
set epl = (`grep :EPL00$iatom $file.scf2$updn`)
set eph = (`grep :EPH00$iatom $file.scf2$updn`)
  else if ($iatom <= 99) then
set epl = (`grep :EPL0$iatom $file.scf2$updn`)
set eph = (`grep :EPH0$iatom $file.scf2$updn`)
  else
set epl = (`grep :EPL$iatom $file.scf2$updn`)
set eph = (`grep :EPH$iatom $file.scf2$updn`)
  endif
if($#epl == 1 ) then
  echo " EPLxx not found in     $file.scf2$updn"
  exit(3)
endif 
if($#eph == 1 ) then
  echo " EPHxx not found in     $file.scf2$updn"
  exit(3)
endif 

#E-parameter will be set when the respective charge is larger than qlimit
unset eslo
unset eplo
unset edlo
unset eflo

if(`echo "if ($epl[2] > $qlimit)1;if ($epl[2] <= $qlimit)0"|bc`) set eslo=$epl[3]
if(`echo "if ($epl[4] > $qlimit)1;if ($epl[4] <= $qlimit)0"|bc`) set eplo=$epl[5]
if(`echo "if ($epl[6] > $qlimit)1;if ($epl[6] <= $qlimit)0"|bc`) set edlo=$epl[7]
if(`echo "if ($epl[8] > $qlimit)1;if ($epl[8] <= $qlimit)0"|bc`) set eflo=$epl[9]
unset es
unset ep
unset ed
unset ef
if(`echo "if ($eph[2] > $qlimit)1;if ($eph[2] <= $qlimit)0"|bc`) set es=$eph[3]
if(`echo "if ($eph[4] > $qlimit)1;if ($eph[4] <= $qlimit)0"|bc`) set ep=$eph[5]
if(`echo "if ($eph[6] > $qlimit)1;if ($eph[6] <= $qlimit)0"|bc`) set ed=$eph[7]
if(`echo "if ($eph[8] > $qlimit)1;if ($eph[8] <= $qlimit)0"|bc`) set ef=$eph[9]

# set the high-E parameter if low-E is already defined
if( $?eslo ) set es=$eph[3]
if( $?eplo ) set ep=$eph[5]
if( $?edlo ) set ed=$eph[7]
if( $?eflo ) set ef=$eph[9]

set n=0
if ($?eslo) @ n ++
if ($?eplo) @ n ++
if ($?edlo) @ n ++
if ($?eflo) @ n ++
if ($?es) then
 if ($es == '10.0000') set es=$eferm
 @ n ++
endif
if ($?ep) then
 if ($ep == '10.0000') set ep=$eferm
 @ n ++
endif
if ($?ed) then
 if ($ed == '10.0000') set ed=$eferm
 @ n ++
endif
if ($?ef) then
 if ($ef == '10.0000') set ef=$eferm
 @ n ++
endif
#write   global energy line for each atom
echo " $eferm   $n   $apw0      global e-param with N other choices, napw ">> $file.in1"$cmplx"new
#write l-specific lines
if ($?es) then
  set de=0.000
  if (  `echo "$es < $etest" | bc` == 1 ) set de=0.002
  echo $es $de $apws | awk '{printf (" 0   % 5.3f     %5.3f CONT%2.0f \n",$1,$2,$3)}' >>$file.in1"$cmplx"new
endif
if ($?eslo) then  
  set de=0.000
  if (  `echo "$eslo < $etest" | bc` == 1 ) set de=0.002
  echo $eslo $de $apws | awk '{printf (" 0   % 5.3f     %5.3f CONT%2.0f \n",$1,$2,$3)}' >>$file.in1"$cmplx"new
endif
if ($?ep) then
  set de=0.000
  if (  `echo "$ep < $etest" | bc` == 1 ) set de=0.002
  echo $ep $de $apwp | awk '{printf (" 1   % 5.3f     %5.3f CONT%2.0f \n",$1,$2,$3)}' >>$file.in1"$cmplx"new
endif
if ($?eplo) then
  set de=0.000
  if (  `echo "$eplo < $etest" | bc` == 1 ) set de=0.002
#  if (  `echo "$eplo < $etest" | bc` == 1 ) echo es set with $de
  echo $eplo $de $apwp | awk '{printf (" 1   % 5.3f     %5.3f CONT%2.0f \n",$1,$2,$3)}' >>$file.in1"$cmplx"new
endif
if ($?ed) then
  set de=0.000
  if (  `echo "$ed < $etest" | bc` == 1 ) set de=0.002
  echo $ed $de $apwd | awk '{printf (" 2   % 5.3f     %5.3f CONT%2.0f \n",$1,$2,$3)}' >>$file.in1"$cmplx"new
endif
if ($?edlo) then
  set de=0.000
  if (  `echo "$edlo < $etest" | bc` == 1 ) set de=0.002
  echo $edlo $de $apwd | awk '{printf (" 2   % 5.3f     %5.3f CONT%2.0f \n",$1,$2,$3)}' >>$file.in1"$cmplx"new
endif
if ($?ef) then
  set de=0.000
  if (  `echo "$ef < $etest" | bc` == 1 ) set de=0.002
  echo $ef $de $apwf | awk '{printf (" 3   % 5.3f     %5.3f CONT%2.0f \n",$1,$2,$3)}' >>$file.in1"$cmplx"new
endif
if ($?eflo) then
  set de=0.000
  if (  `echo "$eflo < $etest" | bc` == 1 ) set de=0.002
  echo $eflo $de $apwf | awk '{printf (" 3   % 5.3f     %5.3f CONT%2.0f \n",$1,$2,$3)}' >>$file.in1"$cmplx"new
endif

   @ iatom ++
end

# @ lines ++
#write rest of case.in1 file
set nn=`wc $file.in1"$cmplx"`
@ n = $nn[1] - $lines
tail -$n $file.in1"$cmplx" >>$file.in1"$cmplx"new

exit(0)


help:

cat <<EOF 
write_in1_lapw [-up/dn -c -ql 0.05 ]
write_in1_lapw   writes a case.in1new file using the Energy parameters from 
                 case.scf2[up/dn] (:EPLxx/EPHxx) and case.in1[c]
EOF


More information about the Wien mailing list