Convergence problem

[eliseo]

Hi all,

  I am trying to perform some riper calculations of the unit cell of a metal complex, Fe(II) mononuclear system, with two molecules in the unit cell.
I was trying PBE or SCAN functionals wit  basis = pob-TZVP and jbas  = dhf-TZVP but always I have convergence problems and after 8-9 scf cycles the energy explodes...

              | TOTAL ENERGY         =       -6433.0146972943    |
              | TOTAL ENERGY         =       -6047.3950058257    |
              | TOTAL ENERGY         =       -6377.5263506254    |
              | TOTAL ENERGY         =       -6623.9976990770    |
              | TOTAL ENERGY         =       -6518.7293471434    |
              | TOTAL ENERGY         =       -6574.3035042539    |
              | TOTAL ENERGY         =       -6627.6332493572    |
              | TOTAL ENERGY         =       -6606.1255709050    |
              | TOTAL ENERGY         =       -6294.1403523448    |
              | TOTAL ENERGY         =       -6545.2222578067    |
              | TOTAL ENERGY         =       -6586.3762207693    |
              | TOTAL ENERGY         =       -4334.0223948992    |
              | TOTAL ENERGY         =       -4101.3356070603    |
              | TOTAL ENERGY         =       -4514.5430856829    |
              | TOTAL ENERGY         =       -2667.6707907154    |

  I was trying to increasing the scfdamp parameter but I obtain more or less the same result.
Any advice will be welcome, I add the control file at the end

     best wishes, thanks a lot

            Eliseo




$title
s13 hs scan periodic D
$operating system unix
$symmetry c1
$user-defined bonds    file=coord
$coord    file=coord
$optimize
 internal   off
 redundant  off
 cartesian  on
 global     off
 basis      off
$atoms
fe 1-2                                                                         \
   basis =fe pob-TZVP                                                          \
   jbas  =fe dhf-TZVP
n  3-6,19-20,37-40,49-52,69-70,87-90                                           \
   basis =n pob-TZVP                                                           \
   jbas  =n dhf-TZVP
c  7-8,11-12,15-16,21-24,27-28,31-32,35-36,41-42,45-46,53-54,57-58,61-62,65-66 \
   71-74,77-78,81-82,85-86,91-92,95-96,99-100                                  \
   basis =c pob-TZVP                                                           \
   jbas  =c dhf-TZVP
h  9-10,13-14,17-18,25-26,29-30,33-34,43-44,47-48,55-56,59-60,63-64,67-68,     \
   75-76,79-80,83-84,93-94,97-98,101-102                                       \
   basis =h pob-TZVP                                                           \
   jbas  =h dhf-TZVP
b  103-104,113-114                                                             \
   basis =b pob-TZVP                                                           \
   jbas  =b dhf-TZVP
f  105-112,115-122                                                             \
   basis =f pob-TZVP                                                           \
   jbas  =f dhf-TZVP
$basis    file=basis
$rundimensions
   dim(fock,dens)=1815888
   natoms=122
   nshell=844
   nbf(CAO)=1904
   dim(trafo[SAO<-->AO/CAO])=2096
   rhfshells=2
   nbf(AO)=1808
$uhfmo_alpha   file=alpha
$uhfmo_beta   file=beta
$uhf
$alpha shells
 a       1-332                                  ( 1 )
$beta shells
 a       1-324                                  ( 1 )
$scfiterlimit       230
$thize     0.10000000E-04
$thime        5
$scfdump
$scfintunit
 unit=30       size=0        file=twoint
$scfdiis
$maxcor    500 MiB  per_core
$drvopt
   cartesian  on
   basis      off
   global     off
   hessian    on
   dipole     on
   nuclear polarizability
$interconversion  off
   qconv=1.d-7
   maxiter=25
$coordinateupdate
   dqmax=0.3
   interpolate  on
   statistics    5
$forceupdate
   ahlrichs numgeo=0  mingeo=3 maxgeo=4 modus=<g|dq> dynamic fail=0.3
   threig=0.005  reseig=0.005  thrbig=3.0  scale=1.00  damping=0.0
$forceinit on
   diag=default
$energy    file=energy
$grad    file=gradient
$forceapprox    file=forceapprox
$dft
   functional pbe
   gridsize   5
   radsize    50
$scfconv   6
$scfdamp   start=0.700  step=0.050  min=0.050
$scforbitalshift  closedshell=.05
$jbas    file=auxbasis
$ricore      500
$rij
$periodic 3
$cell
8.4098 8.4731 18.2909 90.  98.575 90.
$kpoints
nkpoints 2 2 1
$riper
lenonly on
sigma 0.01
$optcell
$disp3 bj
$actual step      riper
$end

[antti_karttunen]

Hi,

which TM version you are using? Periodic calculations are much more robust with TM 7.3.

Best,
Antti

[eliseo]

I was using TURBOMOLE V7.3 ( 22142 )

     thanks,

         Eliseo

[uwe]

Hi,

from what I have heard (not just rumors), convergence problems can be a hint that the input structure is not reasonable. Some popular builder for crystal structures export xyz files which are not defined in a way riper expects it. Especially if you build in symmetry and then export to xyz.

Did you check the structure by visual inspection within Turbomole, i.e. by running TmoleX and view the unit cell plus replicates? Does that look reasonable?

You could try to import the structure directly from a cif file (using TmoleX) and check if you see a similar behavior.

If the structure is OK, and if you have tried changing the other typical parameters already, I'd recommend to contact the Turbomole support.

Regards,

Uwe

[eliseo]

Thanks Uwe for the suggestion about the geometry, I had a problem with the unit cell.
Now, the calculation reaches the convergence in 26 cycles but I am trying to optimize
cell parameters and coordinates together and I have the following message:

 Keyword $statpt not found - using default options

        ***************  Stationary point options ******************
        ************************************************************
           Maximum allowed trust radius:           3.000000E-01
           Minimum allowed trust radius:           1.000000E-03
           Initial trust radius:                   1.500000E-01
           GDIIS used if gradient norm <           1.000000E-02
           Number of previous steps for GDIIS:       5
           Hessian update method:                  BFGS
                        *** Convergence criteria ***
           Threshold for energy change:            1.000000E-06
           Threshold for max displacement element: 1.000000E-03
           Threshold for max gradient element :    1.000000E-03
           Threshold for RMS of displacement:      5.000000E-04
           Threshold for RMS of gradient:          5.000000E-04
        ************************************************************


 keyword $statp missing in file <control>

  Keyword $statpt not found - using default options

 <getgrd> : data group $grad  is missing


========================
 internal module stack:
------------------------
    statpt
========================

 error reading energy and gradient in rdcor
 statpt ended abnormally
statpt step ended abnormally
next step = statpt

Any advice about this problem reading energy and gradient??

    best wishes,

             Eliseo

[antti_karttunen]

Hi,

you should remove the option "lenonly on" under $riper. When this option is included, riper skips the calculation of gradients (needed for optimization).

Best,
Antti

[eliseo]

Hi again,

   Thank you Antti. Now, I have a problem with the orbital occupation in the control file I have

alpha shells
 a       1-332                                  ( 1 )
$beta shells
 a       1-324                                  ( 1 )

  all the parameters seem to perform an unrestricted calculation (define creates alpha and beta files....)
but in the job.last

         ------------------ Fractional Occupations ----------------
           UHF calculation with a single Fermi level for both spins
           Fermi level           =             -0.1608418607
           alpha Electrons(/UC)  =            328.0000000000
           beta  Electrons(/UC)  =            328.0000000000
           Current HOMO-LUMO gap =     0.90223E-01
         ----------------------------------------------------------

              +--------------------------------------------------+
              |                SCF iteration    1                |
              +--------------------------------------------------+

          Number of electrons from P*S  =   655.99999999978832

          Numerical integration of the XC term:
          Number of alpha electrons     =   328.00008873660437
          Number of beta  electrons     =   328.00008873660437

and the energies are equal to those of the restricted calculations. Any suggestion??

    best wishes, thanks a lot

                  Eliseo

[antti_karttunen]

Hi,

In your $riper data group, you have the keyword "sigma 0.01" (Fermi smearing). I have not personally checked out the behavior of this keyword, but at least in the molecular case, the use of Fermi smearing typically changes the occupation to the lowest-eneryg occupation. You can force a certain multiplicity with the "desnue" keyword:

The optional keyword desnue can be used within the $riper data group to constrain the number of unpaired electrons. This can be used to force a certain multiplicity in case of an unrestricted calculation, e.g., desnue 0 for singlet and desnue 1 for dublet.

Please see the manual for further information (Chapter 7.2.4)

Best,
Antti

[eliseo]

Hi again,

   using the desnue keyword, I can perform the unrestricted calculations for the periodic systems with two high-spin Fe(II) in the unit cell (desnue . However, in the restricted low-spin calculation the scf convergence is relatively fast but in the unrestricted high-spin I cannot reach the convergence and after some cycles the energy becomes completely crazy:

              | TOTAL ENERGY         =       -6977.7539894147    |
              | TOTAL ENERGY         =       -6948.8891499268    |
              | TOTAL ENERGY         =       -6966.1812162700    |
              | TOTAL ENERGY         =       -7003.9950034516    |
              | TOTAL ENERGY         =       -6996.5229609826    |
              | TOTAL ENERGY         =       -6982.0301324688    |
.....
              | TOTAL ENERGY         =       -6140.3417339576    |
              | TOTAL ENERGY         =       -2265.1918160646    |
              | TOTAL ENERGY         =       -2289.9961494947    |
              | TOTAL ENERGY         =       -4836.2935906069    |
              | TOTAL ENERGY         =        2140.9916183059    |
              | TOTAL ENERGY         =        2595.2900097573    |
              | TOTAL ENERGY         =        2758.9662267957    |
              | TOTAL ENERGY         =       -1287.5215532921    |
              | TOTAL ENERGY         =       11508.6752014039    |
              | TOTAL ENERGY         =       15459.9397246313   

I think that I don't have geometry problems, I tried as starting geometry the  one employed
for the low-spin state and also the experimental high-spin structure. I was trying many scfdamp
and scforbitalshift but there are not significant differences in the convergence.

    best wishes, thanks in advance

            Eliseo

[antti_karttunen]

Hi, did you also check the SCF Convergence FAQ at http://www.turbo-forum.com/index.php/topic,195.0.html?

Your original control file had "$scforbitalshift  closedshell=.05". Better to use automatic shifting. Here are some settings that are still rather "mild":

Code: [Select]

$scfdamp start=4.500 step=0.200 min=0.500
$scforbitalshift automatic 0.4
$scfiterlimit 200

By the way, is there some particular reason why you are using "gridsize 5" and "radsize 50" for DFT grids? Typically m-grids like "gridsize m4" offer good accuracy and performance.

In the end it could also be that your system is just rather difficult to converge with GGA-PBE.

Best,
Antti

[mas]

Hi,

did you finally manage to converge SCF?

[JakubV]

Hello,

I think it is Gaussian smearing, the "sigma 0.01".
Sigma stands for variance parameter of the gaussian (with centre at Fermi level) in energy.

Fermi smearing would correspond to thermal smearing parametrized by temperature.
Fermi smearing is inside the dscf or ridft modules possible.

I have problems with SCF convergence with RIPER too (only TM7.1 accessible to me right now, sadly). I am desperate enough to try to combine(!) both fermi and gaussian smearings.
But is it reasonable? Will both be applied? And thus convolution of these smearings? My calculation runs so far, no error reported.

Jakub

In your $riper data group, you have the keyword "sigma 0.01" (Fermi smearing).

[turbomaster]

Can you somehow provide me your input files to check?

Best Wishes

Marek Sierka