Dear Sterling,
Maybe your problem is caused by the very small radius for the core correction
in the pseudopotential input file. This leads to a very hard pseudopotential,
which is difficult to describe with the real-space grid. Increasing the
meshcutoff reduces the error, but the forces for a meshcutoff of 2000 Ry are
still 10 times higher than for the other atoms.
i think a rule of thumb for the radius for the partial core correction is the
radius where the valence and core charge are about equal. I recommend you to
try to increase the radius to something like that. The core radii in your
input file are also quite small, the O pseudopotential on the siesta webpage
uses 1.15 bohr.
For good results you will still have to increase the meshcutoff to 200-300 Ry.
Best regards,
Wolfgang
On Tuesday 03 July 2007 07:32, Sterling Paramore wrote:
> Dear SIESTA users,
>
> I'm very new to SIESTA, and I'm having a few problems generating
> pseudopotentials (PP) for Oxygen for use in simulations of
> orthosilicic acid (SiO4H4). I'm trying to follow the work of Ma,
> Foster, and Nieminen (2005,JCP 122, p144709) carefully, so I'm using
> GGA with the PBE functional (I've also been checking it with RPBE,
> which supposedly gives better agreement with experiment). I generated
> the PPs for O using the atom program with the following input file:
>
> pe Oxygen
> tm2 2.0 # PS flavor, logder R
> n=O c=rp # Symbol, XC flavor,{ |r|s}
> 0.0
> 1 4 # norbs_core, norbs_valence
> 2 0 2.00 0.00 # 2s2
> 2 1 4.00 0.00 # 2p2
> 3 2 0.00 0.00 # 3d0
> 4 3 0.00 0.00 # 4f0
> 0.60 0.60 0.60 0.60 0.00 0.075
>
> With the above values, I get good agreement between the all-electron
> and pseudopotential orbital eigenvalues. I generated PP for Si and H
> using very similar input files as well. The problem I'm having is
> that when I go to run a SIESTA calculation (e.g, verlet MD or CG
> minimization), the forces on the oxygen atoms are 3-4 ORDERS of
> MAGNITUDE larger than the forces on the Si or H atoms. Obviously,
> this is not correct.
>
> The strangest issue I'm seeing is that if I generate the O PP with CA
> instead of RP (or PB), and then run SIESTA using GGA and PBE, the
> forces on the Os are comparable to those on Si and H (more like what I
> would expect). This suggests that there's some sort of instability
> that is being generated in the PP when using PB. However, I have
> tried varying just about everything I can in the PP generation with no
> luck. I have tried changing the cutoff from 0.40-2.00, using the full
> core charge by setting rcore_flag=-1, not using the new core
> correction at all (via pg), and even tried using the spin versions.
>
>
> The SIESTA input file I've been using is given below. If anyone has
> any ideas what I might be doing wrong, I'd really appreciate any help.
>
> Sincerely,
> Sterling Paramore
>
> #
> ---------------------------------------------------------------------------
>-- # FDF for orthosilicic acid
> #
> #
> ---------------------------------------------------------------------------
>--
>
> SystemName Orthosilicic Acid
> SystemLabel SiO4H4
>
> #
> ### Species identification
> #
>
> NumberOfSpecies 3
> NumberOfAtoms 9
>
> %block ChemicalSpeciesLabel
> 1 14 Si
> 2 8 O
> 3 1 H
> %endblock ChemicalSpeciesLabel
>
> #
> ### DFT Definitions
>
> XC.functional GGA
> XC.authors RPBE
>
> #
> ### Basis Definitions
> #
>
> User-Basis .false.
> #PAO.BasisSize SZ
> PAO.BasisSize DZP
> PAO.BasisType split
>
> PAO.EnergyShift 272 meV
>
>
> #
> #Boundary Conditions
> #
>
> LatticeConstant 0. Ang
>
> %block LatticeVectors
> 1.000 0.000 0.000
> 0.000 1.000 0.000
> 0.000 0.000 1.000
> %endblock LatticeVectors
>
>
> AtomicCoordinatesFormat Ang
> %block AtomicCoordinatesAndAtomicSpecies
> -0.96890 0.87640 0.04850 1
> -1.83470 1.93920 1.08230 2
> 0.22500 1.75900 -0.81420 2
> -0.21160 -0.34630 0.98630 2
> -2.05870 0.15050 -1.06140 2
> -1.49050 -0.43130 -1.62660 3
> -0.28460 2.43880 -1.32300 3
> -1.14020 2.36020 1.64900 3
> 0.42990 0.14010 1.56310 3
> %endblock AtomicCoordinatesAndAtomicSpecies
>
> #
> # MD Definitions
> #
>
> WriteForces .true.
>
> MD.TypeOfRun CG
> MD.NumCGsteps 10
