# Hydrogen diatomic molecule : computation of derivatives
# of the energy, to a very high accuracy
# Datasets 1 to 5 : GS computations at slightly different geometries,
#  for finite-difference analysis of forces, including the target
#  geometry (for dataset 3)
# Step 6 : RF calculation
# Note : this also tests the use of istwfk==1 in RF with
#  istwfk/=1 in the GS case.

 ndtset 6

    xred1  -0.047 0 0   0.04690  0 0
    xred2  -0.047 0 0   0.04695  0 0
    xred3  -0.047 0 0   0.047    0 0
    xred4  -0.047 0 0   0.04705  0 0
    xred5  -0.047 0 0   0.04710  0 0

    xred6  -0.047 0 0   0.047    0 0

#Specific for RF
  rfphon6  1
 rfatpol6  2 2
   rfdir6  1 0 0
    nqpt6  1
     qpt6  0.0 0.0 0.0
  getwfk6  3
   nstep6  18
  diemix6  0.35
  diemac6  1.0

#Common data
 acell 12 10 10
 amu 1.008
 diemac 1.0d0   diemix 0.5d0
 ecut 4.5
 getwfk -1

 kptopt 0
 kpt   3*0.0
 natom  2
 nband 1

 nkpt 1
 nline 3
 nsym 4  ntypat  1

 rprim 1 0 0  0 1 0  0 0 1
 symrel  1  0  0   0  1  0   0  0  1
         1  0  0   0  1  0   0  0 -1
         1  0  0   0 -1  0   0  0  1
         1  0  0   0 -1  0   0  0 -1
 tnons 12*0
 nstep 12
 tolvrs 7.0d-20
 typat  2*1
 wtk  1
 znucl  1.0

 pp_dirpath "$ABI_PSPDIR"
 pseudos "PseudosTM_pwteter/1h.pspnc"

#%%<BEGIN TEST_INFO>
#%% [setup]
#%% executable = abinit
#%% [files]
#%% files_to_test = 
#%%  t33.out, tolnlines = 0, tolabs = 0.000e+00, tolrel = 0.000e+00
#%% [paral_info]
#%% max_nprocs = 1
#%% [extra_info]
#%% keywords = NC, DFPT
#%% authors = Unknown
#%% description = 
#%%   H2 molecule in a big box : compute VERY accurately
#%%   the derivatives of the energy, by both symmetric finite-differences and
#%%   direct computation of forces and 2DTE.
#%%   Also test the interplay between istwfk/=1 in the GS calculation
#%%   and istwfk==1 in the RF calculation (istwfk/=1 is not yet-991020-
#%%   allowed for RF, which is a shame)
#%%   1) Computation of the first-order derivative of the total energy
#%%   With delta(xred)=0.0002, one gets delta(etot)/delta(xred)=-3.145846551
#%%   With delta(xred)=0.0001, one gets delta(etot)/delta(xred)=-3.145836932
#%%   The combination of both results, in a higher-order finite difference
#%%   formula gives -3.145833726 . The direct computation of forces
#%%   at the target geometry gives -3.145833725869 . The agreement is perfect,
#%%   taking into account the "limited" number of digits (10) of the 
#%%   finite-difference result.
#%%   2) Computation of the second-order derivative of the total energy
#%%   With delta(xred)=0.0002, one gets delta(dedt)/delta(xred)=188.73875
#%%   With delta(xred)=0.0001, one gets delta(dedt)/delta(xred)=188.73837
#%%   The combination of both results, in a higher-order finite difference
#%%   formula gives 188.73824613 . The direct computation of 2DTE
#%%   at the target geometry gives 188.73824613046 . The agreement at the
#%%   level of 11 digits is also perfect.
#%%<END TEST_INFO>