Skip to content

Add ddm_tol as convergence threshold for density matrix #4

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Open
JaafarMehrez wants to merge 1 commit into
base: master
Choose a base branch
from
Open

Add ddm_tol as convergence threshold for density matrix #4

Changes from all commits
Commits
Show all changes
1 commit
Select commit
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
21 changes: 14 additions & 7 deletions pyscf/scf/hf.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -47,7 +47,7 @@
if sys.version_info >= (3,):
unicode = str

def kernel(mf, conv_tol=1e-10, conv_tol_grad=None,
def kernel(mf, conv_tol=1e-10, conv_tol_grad=None, ddm_tol=1e-10,
dump_chk=True, dm0=None, callback=None, conv_check=True, **kwargs):
'''kernel: the SCF driver.

Expand All @@ -74,6 +74,8 @@ def kernel(mf, conv_tol=1e-10, conv_tol_grad=None,
converge threshold.
conv_tol_grad : float
gradients converge threshold.
ddm_tol : float
convergence threshold for density matrix.
dump_chk : bool
Whether to save SCF intermediate results in the checkpoint file
dm0 : ndarray
Expand Down Expand Up @@ -195,7 +197,7 @@ def kernel(mf, conv_tol=1e-10, conv_tol_grad=None,

if callable(mf.check_convergence):
scf_conv = mf.check_convergence(locals())
elif abs(e_tot-last_hf_e) < conv_tol and norm_gorb < conv_tol_grad:
elif abs(e_tot-last_hf_e) < conv_tol and norm_gorb < conv_tol_grad and norm_ddm < ddm_tol:
scf_conv = True

if dump_chk:
Expand Down Expand Up @@ -225,10 +227,11 @@ def kernel(mf, conv_tol=1e-10, conv_tol_grad=None,
norm_ddm = numpy.linalg.norm(dm-dm_last)

conv_tol = conv_tol * 10
ddm_tol = ddm_tol * 10
conv_tol_grad = conv_tol_grad * 3
if callable(mf.check_convergence):
scf_conv = mf.check_convergence(locals())
elif abs(e_tot-last_hf_e) < conv_tol or norm_gorb < conv_tol_grad:
elif abs(e_tot-last_hf_e) < conv_tol or norm_gorb < conv_tol_grad or norm_ddm < ddm_tol:
Copy link
Member

@ndattani ndattani Sep 25, 2024

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I see that we are using norm_ddm which is numpy.linalg.norm(ddm - ddm_last), which is the Frobenius norm of "ddm - ddm_last". This is different from CFOUR which prints the largest density difference.

Copy link
Member Author

@JaafarMehrez JaafarMehrez Jan 26, 2025

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I suppose that the (largest density difference) printed out in CFOUR is the largest element-wise change in the density matrix between iterations. Is that what you are referring to? In this case we could define 'max_ddm = numpy.max(numpy.abs(ddm))' and add that to the convergence criteria instead of norm_ddm.

Copy link
Member Author

@JaafarMehrez JaafarMehrez Jan 26, 2025

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I should mention that I have previously tried using the RMS of the density difference with 'rms_ddm = numpy.sqrt(numpy.mean(numpy.square(ddm)))'.

Compared to 'max_ddm' and 'norm_ddm', I was able to achieve tighter convergence using rms_ddm.

Copy link
Member

@ndattani ndattani Jan 26, 2025

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

@JaafarMehrez yes, when I saw "largest density difference" printed, I always assumed that it was the largest element-wise change in the density matrix between the current iteration and the previous one.

I don't know what you mean by "tighter convergence using rms_ddm". How are you measuring how tight the convergence is? Do you mean that rms_ddm converged to a smaller number than max_ddm did when you tried max_ddm? If max_ddm is 10^-16, then all the other ddm elements will be smaller than 10^-16, which means that rms_ddm will also be small.

The main suggestion I had for you here, was to ensure that PySCF can converge SCF calculations in the same way that more established software such as CFOUR can do it. I recall that you had mentioned on Discord that we can already use the same criterion as MRCC uses.

I have found that in CFOUR, the "largest density difference" cannot be converged as much when we use bigger and bigger basis sets, and this makes it very hard to do the types of calculations that I like to do. However, "largest density difference" seems to be an extremely good criterion because if the "largest density difference" is small, then all other elements of ddm must also be small, meaning that rms and other measures are also likely to be small. Often it took me thousands or tens of thousands of iterations to get the "largest density difference" to be as small as I wanted it to be, but when you do get "largest density difference" converged to a small value, I think you can be very sure that the SCF is very well-converged.

I think a combination of these convergence criteria will be the best thing to use in the long run. Maybe we can try to converge the "largest density difference" until a certain number of iterations is reached, then switch to another traditional convergence criterion if "largest density difference" is not converging after a certain point. These thresholds can be defined by the user, but can also be set to decent default values by the programmer.

Copy link
Member

@ndattani ndattani Jan 26, 2025

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Unfortunately there are now "conflicts" that need to be resolved in order for this pull request to get merged into its upstream branch. Since 4 months have passed since this pull request was made, it is not too surprising that scf/hf.py has changed.

scf_conv = True
logger.info(mf, 'Extra cycle E= %.15g delta_E= %4.3g |g|= %4.3g |ddm|= %4.3g',
e_tot, e_tot-last_hf_e, norm_gorb, norm_ddm)
Expand Down Expand Up @@ -1295,7 +1298,7 @@ def as_scanner(mf):

The solver will automatically use the results of last calculation as the
initial guess of the new calculation. All parameters assigned in the
SCF object (DIIS, conv_tol, max_memory etc) are automatically applied in
SCF object (DIIS, conv_tol, ddm_tol, max_memory etc) are automatically applied in
the solver.

Note scanner has side effects. It may change many underlying objects
Expand Down Expand Up @@ -1372,6 +1375,8 @@ class SCF(lib.StreamObject):
converge threshold. Default is 1e-9
conv_tol_grad : float
gradients converge threshold. Default is sqrt(conv_tol)
ddm_tol : float
convergence threshold for density matrix. Default is 1e-10
max_cycle : int
max number of iterations. If max_cycle <= 0, SCF iteration will
be skiped and the kernel function will compute only the total
Expand Down Expand Up @@ -1436,6 +1441,7 @@ class SCF(lib.StreamObject):
'''
conv_tol = getattr(__config__, 'scf_hf_SCF_conv_tol', 1e-9)
conv_tol_grad = getattr(__config__, 'scf_hf_SCF_conv_tol_grad', None)
ddm_tol = getattr(__config__, 'scf_hf_SCF_ddm_tol', 1e-10)
max_cycle = getattr(__config__, 'scf_hf_SCF_max_cycle', 50)
init_guess = getattr(__config__, 'scf_hf_SCF_init_guess', 'minao')

Expand Down Expand Up @@ -1487,7 +1493,7 @@ def __init__(self, mol):
self.opt = None
self._eri = None # Note: self._eri requires large amount of memory

keys = set(('conv_tol', 'conv_tol_grad', 'max_cycle', 'init_guess',
keys = set(('conv_tol', 'conv_tol_grad', 'ddm_tol', 'max_cycle', 'init_guess',
'DIIS', 'diis', 'diis_space', 'diis_start_cycle',
'diis_file', 'diis_space_rollback', 'damp', 'level_shift',
'direct_scf', 'direct_scf_tol', 'conv_check'))
Expand Down Expand Up @@ -1522,6 +1528,7 @@ def dump_flags(self, verbose=None):
log.info('diis_space = %d', self.diis_space)
log.info('SCF conv_tol = %g', self.conv_tol)
log.info('SCF conv_tol_grad = %s', self.conv_tol_grad)
log.info('SCF ddm_tol = %g', self.ddm_tol)
log.info('SCF max_cycles = %d', self.max_cycle)
log.info('direct_scf = %s', self.direct_scf)
if self.direct_scf:
Expand Down Expand Up @@ -1683,14 +1690,14 @@ def scf(self, dm0=None, **kwargs):
self.converged, self.e_tot, \
self.mo_energy, self.mo_coeff, self.mo_occ = \
kernel(self, self.conv_tol, self.conv_tol_grad,
dm0=dm0, callback=self.callback,
self.ddm_tol, dm0=dm0, callback=self.callback,
conv_check=self.conv_check, **kwargs)
else:
# Avoid to update SCF orbitals in the non-SCF initialization
# (issue #495). But run regular SCF for initial guess if SCF was
# not initialized.
self.e_tot = kernel(self, self.conv_tol, self.conv_tol_grad,
dm0=dm0, callback=self.callback,
self.ddm_tol, dm0=dm0, callback=self.callback,
conv_check=self.conv_check, **kwargs)[1]

logger.timer(self, 'SCF', *cput0)
Expand Down