#PbTiO3 in tetragonal P4mm  structure
#Response function calculation for:
#    * rigid-atom elastic tensor
#    * rigid-atom piezoelectric tensor
#    * interatomic force constants at gamma
#    * Born effective charges

   ndtset   3

# Set 1 : Initial self-consistent run

    iscf1   5
  kptopt1   1
  tolvrs1   1.0d-12  #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
   rfdir2   1 1 1    #full set of directions needed
  tolwfr2   1.0d-12  #only wf convergence can be monitored here

# Set 3 : response-function calculations for all needed perturbations

  getddk3  -1
  getwfk3  -2
    iscf3   5
  kptopt3   2        #use time-reversal symmetry only for k points
    nqpt3   1
     qpt3   0 0 0
  rfphon3   1        #do atomic displacement perturbation
 rfelfd3   3
 rfatpol3   1 5      #do for all atoms
  rfstrs3   3        #do strain perturbation
   rfdir3   1 1 1    #the full set of directions is needed
  tolvrs3   1.0d-12  #need reasonable convergence of 1st-order quantities

#Common input data


# acell  COPY RELAXED RESULT FROM PREVIOUS CALCULATION
# Here is a default value, for automatic testing : suppress it and fill the previous line
  acell  2*7.376 7.8544 



  rprim  1.0  0.0  0.0   #hexagonal primitive vectors must be
         0.0  1.0  0.0   #specified with high accuracy to be
         0.0  0.0  1.0   #sure that the symmetry is recognized
                                 #and preserved in the optimization
                                 #process

#Definition of the atom types and atoms
 ntypat   3 
  znucl   82 81 8
  natom   5
  typat   1 2 3 3 3

#Starting approximation for atomic positions in REDUCED coordinates
#based on ideal tetrahedral bond angles

# xred  COPY RELAXED RESULT FROM PREVIOUS CALCULATION
# Here is a set of default values, for automatic testing : suppress it and fill the previous line
  xred  0.0 0.0 0.1206
        0.5 0.5 0.5773
        0.5 0.5 -0.0077
        0.5 0.0 0.5
        0.0 0.5 0.5
#Gives the number of bands, explicitely (do not take the default)
  nband   13              # 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   20.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   8 8 8          # 4x4x4 Monkhorst-Pack grid
nshiftk   1              # Use one copy of grid only (default)
 shiftk   0.5 0.5 0.5    # 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
   iscf   5              # Use conjugate-gradient SCF cycle
  nstep   200             # Maxiumum number of SCF iterations