[Octopus-users] H2O: optimization of geometry -> plateaus and wiggles

[Секирин И.В.]

     Dear colleagues,
     Could you please explain me, what I am doing incorrectly.
I tried to reproduce vibrational spectrum of the H2O molecule.
First of all, we try to optimize the geometry (see input file below;
I used version 2.1.0). Unfortunately, we obtain the message 
that optimal geometry is not achieved. For finding the reason, 
we try to look at E(r) along the line that the program tries 
to move along. In geom_opt.F90, we insert:
      call geom_calc_point(x, f, df)
      do iter=-10,10
        x1 = x - (iter/5)*g_opt%step * df; call geom_calc_point(x1, f1, df1)
        write(*,'(''geom_opt.248: iter,f1= '',i5,1x,f14.8)')iter,f1
      enddo
      stop
     After re-translation and re-running, we obtain the following:
- for ConvAbsDens = 1e-5, E(r) is noise-like with dispersion ~1e-5; 
-   for ConvAbsDens = 1e-7, E(r) is parabola-like with artifacts (plateaus
  and wiggles) whose amplitude is ~1e-7;
- for ConvAbsDens = 1e-8, sometimes SCF is not aschieved (naturally - with
    two digits before decimal point and ~10 meaningfull digits, ConvAbsDens
    is comparable with machine error).
     E(r) for ConvAbsDens = 1e-5 and 1e-7 is presented below 
(can be wieved using gnuplot).
     Unfortunately, it seems natural that with such E(r) finding minimum
is problematic. 
     However, one could expect that for some geometry close to the minimum
d2E/dr2 should be close to those for minimum, and the calculated 
vibrational spectrum shoud be qualitatively correct. Indeed, small
variations of the obtained geometry lead to small changes in the
calculated frequencies (in both cases - in the fourth meaningfull digit), 
but in all cases, unfortunately, we even do not obtain six zero frequencies 
for translations and rotations - the calulated ones are uniformly 
distributed in the 1000..4500 cm-1 range. I.e., even for optimal geometry,
obtaining vibrational frequencies seems to remain a problem.
     With best regards     Igor.
----- inp -----\
SystemName = 'H2O'

Units = "eVA"

%Species
 'O' | 16.0 | spec_ps_tm2 | 1 | 1 | 1
 'H' | 1.0  | spec_ps_tm2 | 1 | 1 | 1
%

rOH=1.09
# %Coordinates
#  'O' | 0.0 | 0.0 | 0.0 | no
#  'H' | rOH/2 | -rOH*sqrt(3)/2 | 0.0 | yes
#  'H' | rOH/2 |  rOH*sqrt(3)/2 | 0.0 | yes
# %
%Coordinates
 'O' | 0.000760 |   0.000000 |  0.000000 | yes
 'H' | 0.580104 |  -0.774500 |  0.000000 | yes
 'H' | 0.580104 |   0.774500 |  0.000000 | yes
%

Radius = 6.0
Spacing = .2
TypeOfMixing = 2
ElectronicTemperature = 0.1
MaximumIter = 50

# --------

# DebugLevel = 1
FlushMessages = yes

# --------

# CalculationMode = gs
# ExtraStates = 1

# --------

# Before CalculationMode = td one should perform CalculationMode = gs !
# One can TDOutput = multipoles + geometry , but it is so by default.
# Results see td.general/multipoles .

# CalculationMode = td
TDDeltaStrength = 0.05
%TDPolarization
 1 |  1 |  1
 2 | -1 | -1
 0 |  1 | -1
%
TDPolarizationDirection = 1
TDMaximumIter = 20
OutputEvery = 50

# --------

CalculationMode = geom
GOStep = 10.0
GOTolerance = 0.0001
ConvAbsDens = 1e-07

# CalculationMode = phonons

----- inp -----/
----- E(r) for ConvAbsDens = 1e-05 -----\
-20 -17.17976137
-19 -17.17976139
-18 -17.17975998
-17 -17.17976001
-16 -17.17976000
-15 -17.17976002
-14 -17.17975007
-13 -17.17977807
-12 -17.17975737
-11 -17.17976952
-10 -17.17974620
-9 -17.17975485
-8 -17.17977324
-7 -17.17974187
-6 -17.17976747
-5 -17.17975290
-4 -17.17976648
-3 -17.17975250
-2 -17.17978158
-1 -17.17975734
0 -17.17977135
1 -17.17974573
2 -17.17976710
3 -17.17975415
4 -17.17977755
5 -17.17976141
6 -17.17976223
7 -17.17976222
8 -17.17976222
9 -17.17976221
10 -17.17975133
11 -17.17977491
12 -17.17975759
13 -17.17976923
14 -17.17974735
15 -17.17974891
16 -17.17977688
17 -17.17975605
18 -17.17976876
19 -17.17974426
20 -17.17975112
----- E(r) for ConvAbsDens = 1e-05 -----/
----- E(r) for ConvAbsDens = 1e-07 -----\
-20 -17.17975831
-19 -17.17975999
-18 -17.17976010
-17 -17.17976010
-16 -17.17976011
-15 -17.17976011
-14 -17.17976175
-13 -17.17976142
-12 -17.17976143
-11 -17.17976142
-10 -17.17976143
-9 -17.17976197
-8 -17.17976213
-7 -17.17976212
-6 -17.17976214
-5 -17.17976213
-4 -17.17976250
-3 -17.17976235
-2 -17.17976240
-1 -17.17976236
0 -17.17976240
1 -17.17976238
2 -17.17976241
3 -17.17976238
4 -17.17976241
5 -17.17976245
6 -17.17976218
7 -17.17976219
8 -17.17976218
9 -17.17976220
10 -17.17976154
11 -17.17976142
12 -17.17976146
13 -17.17976143
14 -17.17976146
15 -17.17976001
16 -17.17976020
17 -17.17976025
18 -17.17976022
19 -17.17976025
20 -17.17975852
----- E(r) for ConvAbsDens = 1e-07 -----/