TURBOMOLE Users Forum
TURBOMOLE Modules => Riper => Topic started by: eliseo on February 18, 2019, 12:06:11 PM
-
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
-
Hi,
which TM version you are using? Periodic calculations are much more robust with TM 7.3.
Best,
Antti
-
I was using TURBOMOLE V7.3 ( 22142 )
thanks,
Eliseo
-
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
-
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
-
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
-
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
-
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
-
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 8). 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
-
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":
$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
-
Hi,
did you finally manage to converge SCF?
-
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).
-
Can you somehow provide me your input files to check?
Best Wishes
Marek Sierka