GAPW and GTH
Harald Forbert
When calculating HCl with GTH pseudos, I get weird energies with
GAPW (positive total energy: +11.6) whereas with GPW
it is reasonably bound at -15.5. Now if I turn on FULL_GAPW
it again gives sort of the same energy as GPW, but in both GAPW
cases the conserved quantity is not conserved - even in a few (5)
steps MD I get a change of 1e-3 whereas the corresponding
GPW calculation the conserved quantity only changes by 1.0e-7.
I did not do a longer MD with FULL_GAPW, but a similar
calculation with "normal" GAPW also blew up, so it's not just
a miscalculated conserved quantity.
Both GPW+GTH and GAPW with all electrons seem to work fine.
I used the GTH_BASIS_SETS and GTH_POTENTIALS from
the cvs tree and I pasted the input file below.
Am I doing something wrong, or is GAPW not fully working?
The same sort of setup works perfectly fine for a single water
molecule.....
&GLOBAL
PROGRAM_NAME CP2K
PROJECT gapw
RUN_TYPE MD
PRINT_LEVEL LOW
&END GLOBAL
&MOTION
&MD
ENSEMBLE NVE
STEPS 5
TIMESTEP 0.2
TEMPERATURE 50.0
&END MD
&END MOTION
&FORCE_EVAL
METHOD QS
&DFT
BASIS_SET_FILE_NAME GTH_BASIS_SETS
POTENTIAL_FILE_NAME GTH_POTENTIALS
&QS
EPS_DEFAULT 1.0E-14
EXTRAPOLATION PS
EXTRAPOLATION_ORDER 3
METHOD GAPW
# FULL_GAPW T
MAP_CONSISTENT T
&END QS
&SCF
MAX_SCF 100
EPS_SCF 1.0E-7
CHOLESKY F
SCF_GUESS ATOMIC
&END SCF
&MGRID
CUTOFF 400
&END MGRID
&XC
&XC_FUNCTIONAL BLYP
&END XC_FUNCTIONAL
&END XC
&END DFT
&SUBSYS
&CELL
ABC 8.0 8.0 8.0
UNIT ANGSTROM
&END CELL
&COORD
Cl 0.9564675869 -0.4777215658 0.1818080168
H -0.1657168075 0.1752674377 -0.0159958272
&END COORD
&VELOCITY
1.57689068E-05 2.13541290E-05 1.89890498E-06
-2.92031628E-04 6.55172505E-05 1.82242571E-04
&END VELOCITY
&KIND Cl
BASIS_SET TZV2P-GTH
MASS 3.4968852720999998E+01
POTENTIAL GTH-BLYP-q7
&END KIND
&KIND H
BASIS_SET TZV2P-GTH
MASS 1.0079469999999999E+00
POTENTIAL GTH-BLYP-q1
&END KIND
&END SUBSYS
&END FORCE_EVAL
Juerg Hutter
there is an unsolved problem with GAPW: depending on the
basis sets involved, results are very dependent on the
EPS_FIT variable. In you case you can set EPS_FIT to
10^-2 to get reasonable results. Also MD looks much better.
The problem is related to the split of the primitive functions
into hard and soft sets. Changing EPS_FIT moves primitives
from one set to the other. In prinicple this should affect
the total energy only slightly as we trade one source of
errors with another.
I will pu this problem very high on my "to be investigated"
list.
best regards
Juerg
----------------------------------------------------------
Juerg Hutter Phone : ++41 44 635 4491
Physical Chemistry Institute FAX : ++41 44 635 6838
University of Zurich E-mail: hut...@pci.uzh.ch
Winterthurerstrasse 190
CH-8057 Zurich, Switzerland
----------------------------------------------------------
Harald Forbert
work really well. Thanks a lot. :)
Though while testing I noticed something 'odd':
If I use H2O instead of HCl (and essentially the same input as above),
the results also look good, the grid effect is quite a bit smaller
than the
default EPS_FIT, MD seems stable, but if I just do a energy_force
calculation of the following configuration:
O 0.0 0.0 0.0
H 0.0 0.59124 -0.77291
H 0.0 0.59124 0.77291
I get a strange force on the oxygen in z direction:
ENERGY| Total FORCE_EVAL ( QS ) energy (a.u.):
-17.211597494294843
FORCES| (a.u.)
O 0.0000000127 -0.0024393452 0.0001621868
H -0.0000000043 0.0012914592 0.0010373482
H 0.0000000041 0.0011673908 -0.0011995375
Now, the force isn't that large, but if I simply switch axis (from z
to x) to:
O 0.0 0.0 0.0
H -0.77291 0.59124 0.0
H 0.77291 0.59124 0.0
I get basically the same energy, and the force on oxygen looks much
better as
one would expect from symmetry:
ENERGY| Total FORCE_EVAL ( QS ) energy (a.u.):
-17.211597494294828
FORCES| (a.u.)
O -0.0000000021 -0.0024393452 0.0000000127
H 0.0011184426 0.0012294258 -0.0000000043
H -0.0011184430 0.0012294243 0.0000000041
How come the different axis give such different results? Is it really
just the ordering
of the loops over axis in the code? I wouldn't have thought the impact
quite that
large from -2.1e9 to 1.6e-4....
Or is this an unlucky extreme example (due to symmetry)?
BTW the only force term that is really different in both calculations
is the vhxc_atom part and
for default EPS_FIT the difference is far less pronounced.
-H
Matthias Krack
On the basis of my experience an EPS_FIT value of 0.01 is by far too large
to obtain reasonable forces even if energies might look ok. Basically,
this is a (dirty) workaround for a long present bug in the GAPW
implementation which pops up whenever no hard basis is present, ie. the
full atomic basis is equivalent to the soft basis as already outlined by
Juerg. The O basis is much harder than the Cl basis set and thus you could
try smaller EPS_FIT values for water (even smaller than the default
EPS_FIT value). Moreover, for GAPW calculation the atomic KIND parameters
RADIAL_GRID, LEBEDEV_GRID, MAX_RAD_LOCAL might have an impact on the
accuracy (rotational invariance) of your results.
Matthias
Dr. Matthias Krack
Computational Science
Department of Chemistry and Applied Biosciences
ETH Zurich
USI-Campus, via G. Buffi 13
6900 Lugano
Switzerland
Phone: +41 (0)58 666 48 05 (direct)
+41 (0)58 666 48 00 (secretary)
Fax: +41 (0)58 666 48 17
Email: kr...@phys.chem.ethz.ch
URL: http://www.phys.chem.ethz.ch