Hi, everyone!
I’m interested in knowing why an SCF calculation for my system is slow to converge when the system itself was optimized using a high level method and a large basis set. The system in question is a molecule that is interacting with a Zn2+ ion.
I hope it’s not against the forum policy to post inputs and snippet results as I think I should be detailed in explaining what’s happening. Kindly find the input for the optimization calculation below for clarity:
#! Optimization of dimer b3lyp-d/quad zeta with diffuse and polarization
set_num_threads(8)
molecule {
2 1
--
0 1
S 3.9850041190 1.7475301917 -0.1954783522
O 4.1932869082 1.7659313800 1.5064078660
O 5.3009547141 2.2352305701 -1.1827981411
O 2.4631929821 1.5248049343 -0.9516507767
C -2.8731149350 -0.5684502801 0.5725428690
C -4.0896556319 -0.7180563519 -0.3586099991
C -1.5561262291 -0.8834079246 -0.1582708978
C -5.4059608460 -0.4117836361 0.3785907810
C -0.3437899300 -0.7365903296 0.7780077229
C -6.6256497183 -0.5800970287 -0.5448878418
C 0.9722090029 -1.0539977649 0.0488187326
C -7.9399726807 -0.2746319122 0.1943260327
C 2.1941637808 -0.7656488287 0.9418442604
C -9.1590633540 -0.4526448306 -0.7259386774
C 3.4981447121 -1.0642555071 0.2211539598
C -10.4684106102 -0.1445804956 0.0160550780
C 4.3829091816 -0.0564081829 -0.2909681019
C -11.6831647552 -0.3234689137 -0.8934457025
C 3.8383268419 -2.4245785135 0.0719660467
C 5.5793338877 -0.4875799161 -0.9163926008
C 5.0223334372 -2.8060904030 -0.5526549257
C 5.8895204514 -1.8398901393 -1.0423655431
H -2.8359296824 0.4734001219 0.9609218947
H -2.9847109482 -1.2654878965 1.4315975645
H -3.9818507060 -0.0210221275 -1.2188154215
H -4.1231030736 -1.7598641657 -0.7462018320
H -1.5914422253 -1.9221244876 -0.5524202252
H -1.4398545575 -0.1812511319 -1.0141272562
H -5.5084814147 -1.1036954677 1.2433492177
H -5.3765864332 0.6327571618 0.7596970344
H -0.3127093799 0.3103104440 1.1530566240
H -0.4503287204 -1.4234531946 1.6460438291
H -6.6552000374 -1.6261068593 -0.9218511659
H -6.5269553796 0.1090262261 -1.4122428002
H 0.9700993694 -2.1139159094 -0.2838332862
H 1.0359640066 -0.4156804009 -0.8550884309
H -7.9147482622 0.7728642057 0.5675955123
H -8.0390314418 -0.9607142444 1.0642515114
H 2.1470896941 0.2909185261 1.2631482411
H 2.1497643917 -1.3749544844 1.8738568266
H -9.0650142439 0.2318032788 -1.5976538877
H -9.1888132260 -1.5010197877 -1.0965965739
H -10.5742122355 -0.8265242532 0.8877822672
H -10.4516028715 0.9042891580 0.3844769508
H -12.6109155763 -0.0927479887 -0.3273989031
H -11.6214770360 0.3633599334 -1.7643290328
H -11.7444021763 -1.3708372128 -1.2585735364
H 3.1857469455 -3.2005841869 0.4545620884
H 6.3001070506 0.2002457244 -1.3266206808
H 5.2716475911 -3.8548336970 -0.6550158206
H 6.8098380439 -2.1378883093 -1.5288485220
Na 5.6120705484 3.2174854309 1.7005216289
--
2 1
Zn 3.4730119981 3.9307810064 -0.1230901354
}
# both molecules are singlets, but one has charge
set optking { opt_coordinates = cartesian }
set {
maxiter 500
basis def2-QZVPPD
scf_type df
guess sad
# basis_guess def2-SV(P)
reference uks
}
optimize('b3lyp-d')
So, I used the optimized geometry from that, which is the following:
geom = psi4.geometry("""
2 1
S 3.691592181079 1.045529796041 0.438604633103
O 3.541481915892 1.209652609225 1.878132736683
O 4.861740319586 1.769502464398 -0.066283240269
O 2.446534573903 1.494263869452 -0.295728929209
C -3.398986350701 -1.049879235358 0.805468402569
C -4.596341024186 -1.248061153617 -0.144232465459
C -2.087015828839 -1.301758920470 0.041406548169
C -5.895069776288 -0.965442567837 0.644580812758
C -0.884532775597 -1.139375625846 0.980380041836
C -7.081543874581 -1.165614119500 -0.331766871472
C 0.429158223201 -1.387889246436 0.235444642647
C -8.370059057679 -0.852733513607 0.476909824474
C 1.635179063533 -1.302380053198 1.207051672359
C -9.549993928040 -1.044368456867 -0.505951446607
C 2.909751057812 -1.622750365424 0.487051942736
C -10.850886123937 -0.720893905999 0.287577006062
C 3.864003145052 -0.680179284982 0.056059202783
C -12.032420204375 -0.900671413588 -0.676049255631
C 3.163540142727 -2.966370815367 0.170382303622
C 4.999775136938 -1.058461614007 -0.651324464979
C 4.292967599201 -3.357693908108 -0.529114478068
C 5.213902038328 -2.403243582741 -0.949206854386
H -3.406335424835 -0.035620269450 1.206338366688
H -3.476229063539 -1.738987441013 1.647064168575
H -4.520085264342 -0.567918738736 -0.992324469303
H -4.605949408194 -2.267719279915 -0.528253340089
H -2.093391145635 -2.309472367355 -0.377296549043
H -2.008682012276 -0.602633212533 -0.793199431898
H -5.982543316629 -1.648322144934 1.488026604897
H -5.889717860831 0.054477355988 1.025758089674
H -0.881049106485 -0.135917887623 1.410919078210
H -0.972726714436 -1.844287111473 1.809532785822
H -7.103181036909 -2.189230715964 -0.700263233090
H -6.994157924607 -0.489542510306 -1.180000159205
H 0.409625804237 -2.374631303126 -0.229150842881
H 0.557904399161 -0.654552330603 -0.562008991552
H -8.345783300562 0.169863134709 0.847781842653
H -8.462577704485 -1.531519476742 1.322141775479
H 1.655058534294 -0.322188933400 1.673582857981
H 1.484321679350 -2.038670560521 1.997661914959
H -9.455829262088 -0.370711319039 -1.355870631065
H -9.581632238293 -2.068302389032 -0.873932222475
H -10.935095370346 -1.394395780349 1.138851076193
H -10.804850405376 0.300722416715 0.661340199985
H -12.945757599076 -0.671618974015 -0.114384693157
H -11.967371255460 -0.219412278388 -1.523035612203
H -12.099612227869 -1.923410888518 -1.043502018142
H 2.448855152313 -3.710220253737 0.494851242312
H 5.715212396337 -0.313912368877 -0.964851020613
H 4.458137574225 -4.402974601337 -0.746450855203
H 6.095999014670 -2.695246174986 -1.499779434233
Na 5.567285215970 2.399943070327 2.193759087524
Zn 2.501951984597 2.818011352358 -1.739867362110
""")
for energy calculations using the following settings:
# establishing list of methods
methods = ['pbe','b97-d','b3lyp','b3lyp-d']
# establishing karslruhe set
basis_sets = ['cc-pvtz','aug-cc-pvtz','def2-TZVPD','def2-TZVPPD']
# define calculation function
def b3lyp(basis, method, num_threads=4):
psi4.set_options({
'maxiter':500,
'reference':'uks',
'basis':basis,
'scf_type':'df',
'guess':'sad',
'basis_guess':'def2-SV(P)'
})
psi4.set_num_threads(num_threads)
ener = psi4.energy(method, molecule = geom)
return ener
ener_dict={}
for xc in methods:
ener_dict[xc]=[]
for bs in basis_sets:
energy=b3lyp(basis = bs, method = xc, num_threads=8)
ener_dict[xc].append(energy)
pass
pass
While the Delta E
gets progressively smaller, it takes very small steps in doing so. For example, this is a snippet from the first round of the for loop (i.e. using pbe/cc-pvtz
):
@DF-UKS iter 363: -3265.06150783944668 1.06693e-08 1.03369e-03 DIIS
@DF-UKS iter 364: -3265.06150784973761 -1.02909e-08 1.03369e-03 DIIS
@DF-UKS iter 365: -3265.06150783980956 9.92804e-09 1.03369e-03 DIIS
@DF-UKS iter 366: -3265.06150784938291 -9.57334e-09 1.03369e-03 DIIS
@DF-UKS iter 367: -3265.06150784014790 9.23501e-09 1.03369e-03 DIIS
@DF-UKS iter 368: -3265.06150784905822 -8.91032e-09 1.03369e-03 DIIS
@DF-UKS iter 369: -3265.06150784046895 8.58927e-09 1.03369e-03 DIIS
@DF-UKS iter 370: -3265.06150784875445 -8.28550e-09 1.03369e-03 DIIS
@DF-UKS iter 371: -3265.06150784076272 7.99173e-09 1.03369e-03 DIIS
@DF-UKS iter 372: -3265.06150784846704 -7.70433e-09 1.03369e-03 DIIS
@DF-UKS iter 373: -3265.06150784103556 7.43148e-09 1.03369e-03 DIIS
@DF-UKS iter 374: -3265.06150784820602 -7.17046e-09 1.03369e-03 DIIS
@DF-UKS iter 375: -3265.06150784128931 6.91671e-09 1.03369e-03 DIIS
@DF-UKS iter 376: -3265.06150784796318 -6.67387e-09 1.03369e-03 DIIS
@DF-UKS iter 377: -3265.06150784152305 6.44013e-09 1.03369e-03 DIIS
@DF-UKS iter 378: -3265.06150784773490 -6.21185e-09 1.03369e-03 DIIS
@DF-UKS iter 379: -3265.06150784174406 5.99084e-09 1.03369e-03 DIIS
@DF-UKS iter 380: -3265.06150784752026 -5.77620e-09 1.03369e-03 DIIS
@DF-UKS iter 381: -3265.06150784194961 5.57066e-09 1.03369e-03 DIIS
@DF-UKS iter 382: -3265.06150784732472 -5.37511e-09 1.03369e-03 DIIS
@DF-UKS iter 383: -3265.06150784214333 5.18139e-09 1.03369e-03 DIIS
@DF-UKS iter 384: -3265.06150784714100 -4.99767e-09 1.03369e-03 DIIS
@DF-UKS iter 385: -3265.06150784231795 4.82305e-09 1.03369e-03 DIIS
@DF-UKS iter 386: -3265.06150784697002 -4.65207e-09 1.03369e-03 DIIS
@DF-UKS iter 387: -3265.06150784247984 4.49018e-09 1.03369e-03 DIIS
@DF-UKS iter 388: -3265.06150784681449 -4.33465e-09 1.03369e-03 DIIS
@DF-UKS iter 389: -3265.06150784263355 4.18095e-09 1.03369e-03 DIIS
@DF-UKS iter 390: -3265.06150784666715 -4.03361e-09 1.03369e-03 DIIS
@DF-UKS iter 391: -3265.06150784277906 3.88809e-09 1.03369e-03 DIIS
@DF-UKS iter 392: -3265.06150784652527 -3.74621e-09 1.03369e-03 DIIS
@DF-UKS iter 393: -3265.06150784291094 3.61433e-09 1.03369e-03 DIIS
@DF-UKS iter 394: -3265.06150784639703 -3.48609e-09 1.03369e-03 DIIS
@DF-UKS iter 395: -3265.06150784303554 3.36149e-09 1.03369e-03 DIIS
@DF-UKS iter 396: -3265.06150784627880 -3.24326e-09 1.03369e-03 DIIS
@DF-UKS iter 397: -3265.06150784314650 3.13230e-09 1.03369e-03 DIIS
@DF-UKS iter 398: -3265.06150784616602 -3.01952e-09 1.03369e-03 DIIS
@DF-UKS iter 399: -3265.06150784325655 2.90947e-09 1.03369e-03 DIIS
@DF-UKS iter 400: -3265.06150784606507 -2.80852e-09 1.03369e-03 DIIS
@DF-UKS iter 401: -3265.06150784335114 2.71393e-09 1.03369e-03 DIIS
@DF-UKS iter 402: -3265.06150784596684 -2.61571e-09 1.03369e-03 DIIS
@DF-UKS iter 403: -3265.06150784344936 2.51748e-09 1.03369e-03 DIIS
@DF-UKS iter 404: -3265.06150784587680 -2.42744e-09 1.03369e-03 DIIS
@DF-UKS iter 405: -3265.06150784353486 2.34195e-09 1.03369e-03 DIIS
Is slow convergence to be expected due to the workflow, resultant geometry, choice of method, basis set, or any other calculation settings? What changes can I make to the calculation so that it converges faster? I would really appreciate input/suggestions from more experienced Psi4 users.
Thank you very much in advance! I hope y’all have a good day.