pyscf / gpu4pyscf

A plugin to use Nvidia GPU in PySCF package
GNU General Public License v3.0
104 stars 18 forks source link

importing gpu4pyscf breaks hessian #130

Open Acetylsalicylsaeure opened 3 months ago

Acetylsalicylsaeure commented 3 months ago

i get an error when trying to calculate the hessian. it works fine with the usual code, but breaks when i import gpu4pyscf.

regular pySCF:

In [1]: from pyscf import gto
   ...: from pyscf.hessian import thermo
   ...: 
   ...: # First compute nuclear Hessian.
   ...: mol = gto.M(
   ...:     atom = '''O    0.   0.       0
   ...:               H    0.   -0.757   0.587
   ...:               H    0.    0.757   0.587''',
   ...:     basis = '631g')
   ...: 
   ...: mf = mol.RKS().run()
   ...: hessian = mf.Hessian().kernel()
converged SCF energy = -75.8179288833364

after importing gpu4pyscf:

In [1]: from gpu4pyscf.dft import rks

In [2]: from pyscf import gto
   ...: from pyscf.hessian import thermo
   ...: 
   ...: # First compute nuclear Hessian.
   ...: mol = gto.M(
   ...:     atom = '''O    0.   0.       0
   ...:               H    0.   -0.757   0.587
   ...:               H    0.    0.757   0.587''',
   ...:     basis = '631g')
   ...: 
   ...: mf = mol.RKS().run()
   ...: hessian = mf.Hessian().kernel()

converged SCF energy = -75.8179288833364
---------------------------------------------------------------------------
AttributeError                            Traceback (most recent call last)
Cell In[2], line 12
      5 mol = gto.M(
      6     atom = '''O    0.   0.       0
      7               H    0.   -0.757   0.587
      8               H    0.    0.757   0.587''',
      9     basis = '631g')
     11 mf = mol.RKS().run()
---> 12 hessian = mf.Hessian().kernel()

File ~/miniconda3/envs/openmm/lib/python3.11/site-packages/gpu4pyscf/hessian/rhf.py:524, in kernel(hessobj, mo_energy, mo_coeff, mo_occ, atmlst)
    521 else:
    522     hessobj.atmlst = atmlst
--> 524 de = hessobj.hess_elec(mo_energy, mo_coeff, mo_occ, atmlst=atmlst)
    525 hessobj.de = de.get() + hessobj.hess_nuc(hessobj.mol, atmlst=atmlst)
    526 mf = hessobj.base

File ~/miniconda3/envs/openmm/lib/python3.11/site-packages/gpu4pyscf/hessian/rhf.py:59, in hess_elec(hessobj, mo_energy, mo_coeff, mo_occ, mo1, mo_e1, h1ao, atmlst, max_memory, verbose)
     57 mo_occ = cupy.asarray(mo_occ)
     58 mo_coeff = cupy.asarray(mo_coeff)
---> 59 de2 = hessobj.partial_hess_elec(mo_energy, mo_coeff, mo_occ, atmlst,
     60                                 max_memory, log)
     61 t1 = log.timer_debug1('hess elec', *t1)
     62 if h1ao is None:

File ~/miniconda3/envs/openmm/lib/python3.11/site-packages/gpu4pyscf/hessian/rks.py:61, in partial_hess_elec(hessobj, mo_energy, mo_coeff, mo_occ, atmlst, max_memory, verbose)
     59 else:
     60     hyb = mf._numint.hybrid_coeff(mf.xc, spin=mol.spin)
---> 61 de2, ej, ek = rhf_hess._partial_hess_ejk(hessobj, mo_energy, mo_coeff, mo_occ,
     62                                          atmlst, max_memory, verbose,
     63                                          abs(hyb) > 1e-10)
     64 de2 += ej - hyb * ek  # (A,B,dR_A,dR_B)
     66 mem_now = lib.current_memory()[0]

File ~/miniconda3/envs/openmm/lib/python3.11/site-packages/gpu4pyscf/hessian/rhf.py:181, in _partial_hess_ejk(hessobj, mo_energy, mo_coeff, mo_occ, atmlst, max_memory, verbose, with_k)
    178 vk1 = vk1.reshape(3,3,nao,nao)
    179 t1 = log.timer_debug1('contracting int2e_ipvip1 for atom %d'%ia, *t1)
--> 181 ej[i0,i0] += contract('xypq,pq->xy', vj1_diag[:,:,p0:p1], dm0[p0:p1])*2
    182 ek[i0,i0] += contract('xypq,pq->xy', vk1_diag[:,:,p0:p1], dm0[p0:p1])
    183 e1[i0,i0] -= contract('xypq,pq->xy', s1aa[:,:,p0:p1], dme0[p0:p1])*2

File ~/miniconda3/envs/openmm/lib/python3.11/site-packages/gpu4pyscf/lib/cutensor.py:151, in contract(pattern, a, b, alpha, beta, out)
    146 def contract(pattern, a, b, alpha=1.0, beta=0.0, out=None):
    147     '''
    148     a wrapper for general tensor contraction
    149     pattern has to be a standard einsum notation
    150     '''
--> 151     return contraction(pattern, a, b, alpha, beta, out=out)

File ~/miniconda3/envs/openmm/lib/python3.11/site-packages/gpu4pyscf/lib/cutensor.py:109, in contraction(pattern, a, b, alpha, beta, out, op_a, op_b, op_c, algo, jit_mode, compute_desc, ws_pref)
    105 alpha = np.asarray(alpha)
    106 beta = np.asarray(beta)
    108 cutensor_backend.contract(cutensor._get_handle().ptr, plan.ptr,
--> 109                          alpha.ctypes.data, a.data.ptr, b.data.ptr,
    110                          beta.ctypes.data, c.data.ptr, out.data.ptr,
    111                          ws.data.ptr, ws_size)
    113 return out

AttributeError: 'memoryview' object has no attribute 'ptr'

i get the same error after just calling to_gpu() method

sunqm commented 3 months ago

Please use gpu4pyscf.dft.RKS(mol) instead of mol.RKS(). The later creates a RKS object implemented in CPU. CPU module and GPU module are not compatible in this case.

There needs a check or a conversion on the ._scf attribute, when initializing hessian, gradients and other modules that rely on mf objects. @wxj6000

wxj6000 commented 3 months ago

Please use gpu4pyscf.dft.RKS(mol) instead of mol.RKS(). The later creates a RKS object implemented in CPU. CPU module and GPU module are not compatible in this case.

There needs a check or a conversion on the ._scf attribute, when initializing hessian, gradients and other modules that rely on mf objects. @wxj6000

This is one of the cases I didn't check the API compatibility carefully. I will work on it this coming week.

@Acetylsalicylsaeure Thank you for your feedback. Please use the solution by Qiming at this point. We are in the process of improving the compatibility of PySCF and GPU4PySCF.

Acetylsalicylsaeure commented 3 months ago

i was more worried about the occurence of pure pyscf code suddenly breaking due to gpu4pyscf being imported using gpu4pyscf.dft.rks.RKS(mol) still leads to the same error

In [1]: import numpy as np
   ...: import pyscf
   ...: from pyscf import lib, df
   ...: from gpu4pyscf.dft import rks
   ...: lib.num_threads(8)
   ...: 
   ...: atom ='''
   ...: O       0.0000000000    -0.0000000000     0.1174000000
   ...: H      -0.7570000000    -0.0000000000    -0.4696000000
   ...: H       0.7570000000     0.0000000000    -0.4696000000
   ...: '''
   ...: 
   ...: mol = pyscf.M(atom=atom, basis='def2-tzvpp', verbose=0)
   ...: 

In [2]: mf = rks.RKS(mol)

In [3]: mf.kernel()
Out[3]: -75.90464116157013

In [4]: hessobj = mf.Hessian()

In [5]: hess = hessobj.kernel()
---------------------------------------------------------------------------
AttributeError                            Traceback (most recent call last)
Cell In[5], line 1
----> 1 hess = hessobj.kernel()

File ~/miniconda3/envs/openmm/lib/python3.11/site-packages/gpu4pyscf/hessian/rhf.py:524, in kernel(hessobj, mo_energy, mo_coeff, mo_occ, atmlst)
    521 else:
    522     hessobj.atmlst = atmlst
--> 524 de = hessobj.hess_elec(mo_energy, mo_coeff, mo_occ, atmlst=atmlst)
    525 hessobj.de = de.get() + hessobj.hess_nuc(hessobj.mol, atmlst=atmlst)
    526 mf = hessobj.base

File ~/miniconda3/envs/openmm/lib/python3.11/site-packages/gpu4pyscf/hessian/rhf.py:59, in hess_elec(hessobj, mo_energy, mo_coeff, mo_occ, mo1, mo_e1, h1ao, atmlst, max_memory, verbose)
     57 mo_occ = cupy.asarray(mo_occ)
     58 mo_coeff = cupy.asarray(mo_coeff)
---> 59 de2 = hessobj.partial_hess_elec(mo_energy, mo_coeff, mo_occ, atmlst,
     60                                 max_memory, log)
     61 t1 = log.timer_debug1('hess elec', *t1)
     62 if h1ao is None:

File ~/miniconda3/envs/openmm/lib/python3.11/site-packages/gpu4pyscf/hessian/rks.py:61, in partial_hess_elec(hessobj, mo_energy, mo_coeff, mo_occ, atmlst, max_memory, verbose)
     59 else:
     60     hyb = mf._numint.hybrid_coeff(mf.xc, spin=mol.spin)
---> 61 de2, ej, ek = rhf_hess._partial_hess_ejk(hessobj, mo_energy, mo_coeff, mo_occ,
     62                                          atmlst, max_memory, verbose,
     63                                          abs(hyb) > 1e-10)
     64 de2 += ej - hyb * ek  # (A,B,dR_A,dR_B)
     66 mem_now = lib.current_memory()[0]

File ~/miniconda3/envs/openmm/lib/python3.11/site-packages/gpu4pyscf/hessian/rhf.py:181, in _partial_hess_ejk(hessobj, mo_energy, mo_coeff, mo_occ, atmlst, max_memory, verbose, with_k)
    178 vk1 = vk1.reshape(3,3,nao,nao)
    179 t1 = log.timer_debug1('contracting int2e_ipvip1 for atom %d'%ia, *t1)
--> 181 ej[i0,i0] += contract('xypq,pq->xy', vj1_diag[:,:,p0:p1], dm0[p0:p1])*2
    182 ek[i0,i0] += contract('xypq,pq->xy', vk1_diag[:,:,p0:p1], dm0[p0:p1])
    183 e1[i0,i0] -= contract('xypq,pq->xy', s1aa[:,:,p0:p1], dme0[p0:p1])*2

File ~/miniconda3/envs/openmm/lib/python3.11/site-packages/gpu4pyscf/lib/cutensor.py:151, in contract(pattern, a, b, alpha, beta, out)
    146 def contract(pattern, a, b, alpha=1.0, beta=0.0, out=None):
    147     '''
    148     a wrapper for general tensor contraction
    149     pattern has to be a standard einsum notation
    150     '''
--> 151     return contraction(pattern, a, b, alpha, beta, out=out)

File ~/miniconda3/envs/openmm/lib/python3.11/site-packages/gpu4pyscf/lib/cutensor.py:109, in contraction(pattern, a, b, alpha, beta, out, op_a, op_b, op_c, algo, jit_mode, compute_desc, ws_pref)
    105 alpha = np.asarray(alpha)
    106 beta = np.asarray(beta)
    108 cutensor_backend.contract(cutensor._get_handle().ptr, plan.ptr,
--> 109                          alpha.ctypes.data, a.data.ptr, b.data.ptr,
    110                          beta.ctypes.data, c.data.ptr, out.data.ptr,
    111                          ws.data.ptr, ws_size)
    113 return out

AttributeError: 'memoryview' object has no attribute 'ptr'

In [6]: 

however, calling .density_fit() on the RKS object fixes the error, is this expected behavior?

Acetylsalicylsaeure commented 3 months ago

Please use gpu4pyscf.dft.RKS(mol) instead of mol.RKS(). The later creates a RKS object implemented in CPU. CPU module and GPU module are not compatible in this case. There needs a check or a conversion on the ._scf attribute, when initializing hessian, gradients and other modules that rely on mf objects. @wxj6000

This is one of the cases I didn't check the API compatibility carefully. I will work on it this coming week.

@Acetylsalicylsaeure Thank you for your feedback. Please use the solution by Qiming at this point. We are in the process of improving the compatibility of PySCF and GPU4PySCF.

ah, i was too slow, thanks! i also experimented around a bit, works well with RHF

wxj6000 commented 3 months ago

@Acetylsalicylsaeure Hessian calculation of regular SCF is not supported on GPU yet. It is only supported for density fitting scheme.

wxj6000 commented 3 months ago

@Acetylsalicylsaeure This is a bug in GPU4PySCF RKS hessian module. PySCF hessian has been redefined. The issue has been fixed in the above PR. Thank you for your feedback.