search

pyscf / pyscf

Python module for quantum chemistry
Apache License 2.0
1.23k stars 572 forks source link

MemoryError in numpy_helper.py #224

Closed zwang123 closed 5 years ago

zwang123 commented 6 years ago

Hi,

I got a MemoryError when performing geometric optimization. The error is at pyscf/lib/numpy_helper.py when creating a large matrix using numpy.empty. Is it possible to check the size of a matrix before creating it and determine whether to store it in memory or on disk?

These are the outputs (including inputs):

from pyscf import gto, dft, scf, cc, geomopt, lib
from berny_aux import * # read_coords
from sys import argv, stderr
import time

log = lib.logger.Logger(stderr, 5)
lib.logger.TIMER_LEVEL = 5

try:
    xyzname = argv[1]
except IndexError:
    xyzname = "h9o4.xyz"

mol = gto.M(
          atom=read_coords(xyzname, format="xyz")
        , charge=1
        , verbose=4
        , basis='augccpvtz')

t0 = time.clock()
mycc = cc.CCSD(scf.RHF(mol))
opt_mol = geomopt.optimize(mycc, assert_convergence=True, 
        verbose=4)

INFO: **** input file end ****

System: uname_result(system='Linux', node='c506-032.stampede2.tacc.utexas.edu', release='3.10.0-693.17.1.el7.x86_64', version='#1 SMP Thu Jan 25 20:13:58 UTC 2018', machine='x86_64', processor='x86_64') Threads 96 Python 3.7.0 (default, Jul 12 2018, 18:15:30) [GCC Intel(R) C++ gcc 6.3 mode] numpy 1.14.5 scipy 1.1.0 Date: Sat Sep 1 12:51:18 2018 PySCF version 1.5.2 PySCF path $$$$$$$$

[ENV] PYSCF_MAX_MEMORY 172800 [ENV] PYSCF_TMPDIR $$$$$$$$ [CONFIG] conf_file None [INPUT] verbose = 4 [INPUT] num. atoms = 13 [INPUT] num. electrons = 40 [INPUT] charge = 1 [INPUT] spin (= nelec alpha-beta = 2S) = 0 [INPUT] symmetry False subgroup None [INPUT] Mole.unit = angstrom [INPUT] 1 O 7.657450000000 16.613500000000 7.094980000000 AA 14.470483312551 31.394964970462 13.407569059267 Bohr [INPUT] 2 H 7.279450000000 17.552000000000 7.193710000000 AA 13.756166837465 33.168472938366 13.594141719545 Bohr [INPUT] 3 H 8.192890000000 16.440000000000 6.126420000000 AA 15.482318268688 31.067097487850 11.577255924058 Bohr [INPUT] 4 H 6.929560000000 15.817500000000 7.003100000000 AA 13.094970563741 29.890742975308 13.233941022942 Bohr [INPUT] 5 O 5.906270000000 14.942500000000 6.978530000000 AA 11.161232717735 28.237232616313 13.187510452061 Bohr [INPUT] 6 H 5.203250000000 14.841700000000 6.310620000000 AA 9.832717457643 28.046748222957 11.925343476203 Bohr [INPUT] 7 H 5.443430000000 14.996100000000 7.872190000000 AA 10.286591878241 28.338521936590 14.876283100540 Bohr [INPUT] 8 O 6.413880000000 18.704900000000 7.619550000000 AA 12.120476595825 35.347138187377 14.398862692430 Bohr [INPUT] 9 H 5.546280000000 18.907000000000 7.187560000000 AA 10.480950210153 35.729051837152 13.582519903879 Bohr [INPUT] 10 H 6.305690000000 18.289100000000 8.555740000000 AA 11.916027126409 34.561390064783 16.168005392986 Bohr [INPUT] 11 O 8.813210000000 16.358200000000 4.744360000000 AA 16.654553178278 30.912517890860 8.965541036341 Bohr [INPUT] 12 H 8.226550000000 16.040600000000 3.959720000000 AA 15.545926450041 30.312340873698 7.482786329963 Bohr [INPUT] 13 H 9.735590000000 16.644900000000 4.466330000000 AA 18.397598761054 31.454302370773 8.440140481929 Bohr

nuclear repulsion = 146.323241552806 number of shells = 133 number of NR pGTOs = 457 number of NR cGTOs = 391 basis = augccpvtz ecp = {} CPU time: 0.98 Create scanner for <class 'pyscf.grad.ccsd.Gradients'> Set <class 'pyscf.cc.ccsd.CCSD'> as a scanner Create scanner for <class 'pyscf.scf.hf.RHF'> 0 Internal coordinates: 0 Number of fragments: 3 0 Number of internal coordinates: 180 0 Number of strong bonds: 10 0 Number of weak bonds: 16 0 Number of strong angles: 9 0 Number of weak angles: 43 0 Number of superweak angles: 15 0 Number of strong dihedrals: 4 0 * Number of weak dihedrals: 83 New geometry (unit Bohr) 1 O 14.470483312551 31.394964970462 13.407569059267 2 H 13.756166837465 33.168472938366 13.594141719545 3 H 15.482318268688 31.067097487850 11.577255924058 4 H 13.094970563741 29.890742975308 13.233941022942 5 O 11.161232717735 28.237232616313 13.187510452061 6 H 9.832717457643 28.046748222957 11.925343476203 7 H 10.286591878241 28.338521936590 14.876283100540 8 O 12.120476595825 35.347138187377 14.398862692430 9 H 10.480950210153 35.729051837152 13.582519903879 10 H 11.916027126409 34.561390064783 16.168005392986 11 O 16.654553178278 30.912517890860 8.965541036341 12 H 15.545926450041 30.312340873698 7.482786329963 13 H 18.397598761054 31.454302370773 8.440140481929

**** <class 'pyscf.scf.hf.as_scanner..SCF_Scanner'> flags **** method = SCF_Scanner-RHF initial guess = minao damping factor = 0 level shift factor = 0 DIIS = <class 'pyscf.scf.diis.CDIIS'> DIIS start cycle = 1 DIIS space = 8 SCF tol = 1e-09 SCF gradient tol = None max. SCF cycles = 50 direct_scf = True direct_scf_tol = 1e-13 chkfile to save SCF result = $$$$$$ max_memory 172800 MB (current use 113 MB) Set gradient conv threshold to 3.16228e-05 init E= -304.226002898833 HOMO = -0.455234142142364 LUMO = -0.0101612366390812 cycle= 1 E= -304.483693139399 delta_E= -0.258 |g|= 0.658 |ddm|= 1.76 HOMO = -0.638153043458442 LUMO = -0.114906123215324 cycle= 2 E= -304.553762286779 delta_E= -0.0701 |g|= 0.295 |ddm|= 0.573 HOMO = -0.700622626139862 LUMO = -0.111276807550624 cycle= 3 E= -304.567231842377 delta_E= -0.0135 |g|= 0.0748 |ddm|= 0.17 HOMO = -0.685539906372402 LUMO = -0.108388755062248 cycle= 4 E= -304.568292301015 delta_E= -0.00106 |g|= 0.0153 |ddm|= 0.0489 HOMO = -0.688361517562838 LUMO = -0.108600531051493 cycle= 5 E= -304.56835313062 delta_E= -6.08e-05 |g|= 0.00534 |ddm|= 0.0109 HOMO = -0.688621305785161 LUMO = -0.108585768275235 cycle= 6 E= -304.568360624741 delta_E= -7.49e-06 |g|= 0.00182 |ddm|= 0.00422 HOMO = -0.688753893859935 LUMO = -0.108570311144226 cycle= 7 E= -304.568361659865 delta_E= -1.04e-06 |g|= 0.000263 |ddm|= 0.00165 HOMO = -0.688812589948163 LUMO = -0.108571161722535 cycle= 8 E= -304.568361700251 delta_E= -4.04e-08 |g|= 7.42e-05 |ddm|= 0.000424 HOMO = -0.688806586644099 LUMO = -0.108570850060705 cycle= 9 E= -304.568361702165 delta_E= -1.91e-09 |g|= 3.15e-05 |ddm|= 8.63e-05 HOMO = -0.688810996214512 LUMO = -0.108570367240245 cycle= 10 E= -304.568361702442 delta_E= -2.77e-10 |g|= 7.19e-06 |ddm|= 3.09e-05 HOMO = -0.688810508365556 LUMO = -0.108570399674339 Extra cycle E= -304.568361702463 delta_E= -2.14e-11 |g|= 3.75e-06 |ddm|= 8.71e-06 converged SCF energy = -304.568361702463

**** <class 'pyscf.cc.ccsd.as_scanner..CCSD_Scanner'> flags **** CC2 = 0 CCSD nocc = 20, nmo = 391 max_cycle = 50 direct = 0 conv_tol = 1e-07 conv_tol_normt = 1e-05 diis_space = 6 diis_start_cycle = 0 diis_start_energy_diff = 1e+09 max_memory 172800 MB (current use 25418 MB) Init t2, MP2 energy = -1.16329686142831 Init E(CCSD) = -1.16329686143262 cycle = 1 E(CCSD) = -1.15258473709442 dE = 0.0107121243 norm(t1,t2) = 0.0601292 cycle = 2 E(CCSD) = -1.17610790290594 dE = -0.0235231658 norm(t1,t2) = 0.0179682 cycle = 3 E(CCSD) = -1.17642009018261 dE = -0.000312187277 norm(t1,t2) = 0.00817064 cycle = 4 E(CCSD) = -1.17867731280743 dE = -0.00225722262 norm(t1,t2) = 0.00299908 cycle = 5 E(CCSD) = -1.17886863597328 dE = -0.000191323166 norm(t1,t2) = 0.000671718 cycle = 6 E(CCSD) = -1.17877295186934 dE = 9.56841039e-05 norm(t1,t2) = 0.000227608 cycle = 7 E(CCSD) = -1.17875252833902 dE = 2.04235303e-05 norm(t1,t2) = 7.72549e-05 cycle = 8 E(CCSD) = -1.17876147629979 dE = -8.94796077e-06 norm(t1,t2) = 3.28041e-05 cycle = 9 E(CCSD) = -1.17876119473506 dE = 2.81564738e-07 norm(t1,t2) = 1.5282e-05 cycle = 10 E(CCSD) = -1.17876090343503 dE = 2.91300023e-07 norm(t1,t2) = 5.49982e-06 cycle = 11 E(CCSD) = -1.17876091636877 dE = -1.29337359e-08 norm(t1,t2) = 2.28e-06 CCSD_Scanner converged E(CCSD_Scanner) = -305.747122618832 E_corr = -1.178760916368768 cycle = 1 norm(lambda1,lambda2) = 0.0163951 cycle = 2 norm(lambda1,lambda2) = 0.00311351 cycle = 3 norm(lambda1,lambda2) = 0.00186498 cycle = 4 norm(lambda1,lambda2) = 0.00068405 cycle = 5 norm(lambda1,lambda2) = 0.00014128 cycle = 6 norm(lambda1,lambda2) = 4.89159e-05 cycle = 7 norm(lambda1,lambda2) = 1.38974e-05 cycle = 8 norm(lambda1,lambda2) = 5.5754e-06 Traceback (most recent call last): File "opt_h9o4.py", line ##, in verbose=4) File "###/pyscf-1.5.2/pyscf/geomopt/init.py", line 19, in optimize return berny_solver.optimize(method, *args, kwargs) File "###/pyscf-1.5.2/pyscf/geomopt/berny_solver.py", line 160, in optimize energy, gradients = solver.send(geom) File "###/pyscf-1.5.2/pyscf/geomopt/berny_solver.py", line 95, in as_berny_solver energy, gradients = g_scanner(mol) File "###/pyscf-1.5.2/pyscf/grad/ccsd.py", line 247, in call mf_grad=mf_grad, kwargs) File "###/pyscf-1.5.2/pyscf/grad/ccsd.py", line 435, in kernel mf_grad, verbose=log) File "###/pyscf-1.5.2/pyscf/grad/ccsd.py", line 64, in kernel d2 = ccsd_rdm._gamma2_outcore(mycc, t1, t2, l1, l2, fdm2, True) File "###/pyscf-1.5.2/pyscf/cc/ccsd_rdm.py", line 163, in _gamma2_outcore gvvvv += lib.einsum('jabc,jd->abcd', jabc, t1) File "###/pyscf-1.5.2/pyscf/lib/numpy_helper.py", line 185, in einsum return dot(At,Bt).reshape(shapeCt, order='A').transpose(new_orderCt) File "###/pyscf-1.5.2/pyscf/lib/numpy_helper.py", line 665, in dot return ddot(a, b, alpha, c, beta) File "###/pyscf-1.5.2/pyscf/lib/numpy_helper.py", line 621, in ddot c = numpy.empty((m,n)) MemoryError

Thank you.

sunqm commented 6 years ago

It's a problem of memory use estimation in CCSD density matrix code. Some memory may not be fully released the moment that the python object was deleted. It causes small overhead in memory usage. This might be the reason of the MemoryError. I occasionally observed this issue when max_memory was set to a value closed to the physical ram size.

The latest fix updated the memory estimation code and a conservative estimation was applied. However, the memory overhead may still be there. A safe solution is to decrease the "max_memory" for the CCSD code a little bit, to leave more spare memory as the buffer for the overhead.

zwang123 commented 6 years ago

I tried a conservative MAX_MEMORY, and got this error message:

Traceback (most recent call last):
  File "###.py", line 64, in <module>
    verbose=4)
  File "###/pyscf/geomopt/__init__.py", line 19, in optimize
    return berny_solver.optimize(method, *args, **kwargs)
  File "###/pyscf/geomopt/berny_solver.py", line 160, in optimize
    energy, gradients = solver.send(geom)
  File "###/pyscf/geomopt/berny_solver.py", line 95, in as_berny_solver
    energy, gradients = g_scanner(mol)
  File "###/pyscf/grad/ccsd.py", line 247, in __call__
    mf_grad=mf_grad, **kwargs)
  File "###/pyscf/grad/ccsd.py", line 435, in kernel
    mf_grad, verbose=log)
  File "###/pyscf/grad/ccsd.py", line 78, in kernel
    _rdm2_mo2ao(mycc, d2, mo_active, fdm2)  # transform the active orbitals
  File "###/pyscf/grad/ccsd.py", line 330, in _rdm2_mo2ao
    fswap.create_dataset('v', (nao_pair,nvir_pair), 'f8', chunks=(nao_pair,blksize))
  File "###/lib/python3.7/site-packages/h5py/_hl/group.py", line 116, in create_dataset
    dsid = dataset.make_new_dset(self, shape, dtype, data, **kwds)
  File "###/lib/python3.7/site-packages/h5py/_hl/dataset.py", line 140, in make_new_dset
    dset_id = h5d.create(parent.id, None, tid, sid, dcpl=dcpl)
  File "h5py/_objects.pyx", line 54, in h5py._objects.with_phil.wrapper
  File "h5py/_objects.pyx", line 55, in h5py._objects.with_phil.wrapper
  File "h5py/h5d.pyx", line 79, in h5py.h5d.create
ValueError: Unable to create dataset (chunk size must be < 4GB)

HDF5 restrains the chunk size to < 4GB. (https://support.hdfgroup.org/HDF5/faq/limits.html)

sunqm commented 6 years ago

The problem has been fixed in pyscf-1.5.3 (and higher).

For older version, a workaround is to reduce the size of max_memory.