NLESC-JCER / QMCTorch

Pytorch Implementation of Real Space Quantum Monte Carlo Simulations of Molecular Systems
https://qmctorch.readthedocs.io/
Apache License 2.0
23 stars 2 forks source link

Inconsistent dimensons #150

Closed scemama closed 8 months ago

scemama commented 10 months ago

Hi! I have tried to follow the documentation for the water molecule, but I used the cc-pvdz basis set (see my bug report you fixed recently). Now, I get an error about inconsistent dimensions in tensors. Note: with this basis, the number of AOs is 25 and the number of MOs is 24, so the coefficient matrix is not square. It is possible that you used the number of MOs instead of the number of AOs somewhere...

In [1]: import numpy as np
In [2]: import matplotlib.pyplot as plt
In [3]: from qmctorch.scf import Molecule
INFO:QMCTorch|  ____    __  ______________             _
INFO:QMCTorch| / __ \  /  |/  / ___/_  __/__  ________/ /  
INFO:QMCTorch|/ /_/ / / /|_/ / /__  / / / _ \/ __/ __/ _ \ 
INFO:QMCTorch|\___\_\/_/  /_/\___/ /_/  \___/_/  \__/_//_/ 
In [5]: from qmctorch.wavefunction import SlaterJastrow
In [6]: from qmctorch.sampler import Metropolis
In [7]: from qmctorch.solver import Solver
In [8]: from qmctorch.utils import plot_walkers_traj
In [10]: mol = Molecule(atom='water.xyz', unit='angs', calculator='pyscf', basis='cc-pvdz', n
    ...: ame='water')
INFO:QMCTorch|
INFO:QMCTorch| SCF Calculation
INFO:QMCTorch|  Running scf  calculation
converged SCF energy = -76.0267987066294
INFO:QMCTorch|  Molecule name       : water
INFO:QMCTorch|  Number of electrons : 10
INFO:QMCTorch|  SCF calculator      : pyscf
INFO:QMCTorch|  Basis set           : cc-pvdz
INFO:QMCTorch|  SCF                 : HF
INFO:QMCTorch|  Number of AOs       : 25
INFO:QMCTorch|  Number of MOs       : 24
INFO:QMCTorch|  SCF Energy          : -76.027 Hartree

In [11]: wf = SlaterJastrow(mol, configs='ground_state')
    ...: 
INFO:QMCTorch|
INFO:QMCTorch| Wave Function
INFO:QMCTorch|  Jastrow factor      : True
INFO:QMCTorch|  Jastrow kernel      : PadeJastrowKernel
INFO:QMCTorch|  Highest MO included : 24
INFO:QMCTorch|  Configurations      : ground_state
INFO:QMCTorch|  Number of confs     : 1
INFO:QMCTorch|  Kinetic energy      : jacobi
INFO:QMCTorch|  Number var  param   : 628
INFO:QMCTorch|  Cuda support        : False

In [12]: sampler = Metropolis(nwalkers=100, nstep=500, step_size=0.25,
    ...:                      nelec=wf.nelec, ndim=wf.ndim,
    ...:                      init=mol.domain('atomic'),
    ...:                      move={'type': 'all-elec', 'proba': 'normal'})
    ...: 
INFO:QMCTorch|
INFO:QMCTorch| Monte-Carlo Sampler
INFO:QMCTorch|  Number of walkers   : 100
INFO:QMCTorch|  Number of steps     : 500
INFO:QMCTorch|  Step size           : 0.25
INFO:QMCTorch|  Thermalization steps: -1
INFO:QMCTorch|  Decorelation steps  : 1
INFO:QMCTorch|  Walkers init pos    : atomic
INFO:QMCTorch|  Move type           : all-elec
INFO:QMCTorch|  Move proba          : normal

In [13]: solver = Solver(wf=wf, sampler=sampler)
    ...: 
INFO:QMCTorch|
INFO:QMCTorch| QMC Solver 
INFO:QMCTorch|  WaveFunction        : SlaterJastrow
INFO:QMCTorch|  Sampler             : Metropolis

In [14]: pos = sampler(wf.pdf)
    ...: pos = pos.reshape(100,10,3).cpu().detach().numpy()
    ...: plt.scatter(pos[:,:,0],pos[:,:,1],s=0.5)
    ...: 
---------------------------------------------------------------------------
RuntimeError                              Traceback (most recent call last)
<ipython-input-14-2c734dc3c3d8> in <module>()
----> 1 pos = sampler(wf.pdf)
      2 pos = pos.reshape(100,10,3).cpu().detach().numpy()
      3 plt.scatter(pos[:,:,0],pos[:,:,1],s=0.5)

/home/scemama/QMCTorch/QMCTorch/qmctorch/sampler/metropolis.py in __call__(self, pdf, pos, with_tqdm)
     93 
     94             self.walkers.initialize(pos=pos)
---> 95             fx = pdf(self.walkers.pos)
     96 
     97             fx[fx == 0] = eps

/home/scemama/QMCTorch/QMCTorch/qmctorch/wavefunction/wf_base.py in pdf(self, pos, return_grad)
    224             return self.gradients(pos, pdf=True)
    225         else:
--> 226             return (self.forward(pos)**2).reshape(-1)
    227 
    228     def get_number_parameters(self):

/home/scemama/QMCTorch/QMCTorch/qmctorch/wavefunction/slater_jastrow.py in forward(self, x, ao)
     92         # atomic orbital
     93         if ao is None:
---> 94             x = self.ao(x)
     95         else:
     96             x = ao

/home/scemama/miniconda3/lib/python3.7/site-packages/torch/nn/modules/module.py in _call_impl(self, *input, **kwargs)
   1192         if not (self._backward_hooks or self._forward_hooks or self._forward_pre_hooks or _global_backward_hooks
   1193                 or _global_forward_hooks or _global_forward_pre_hooks):
-> 1194             return forward_call(*input, **kwargs)
   1195         # Do not call functions when jit is used
   1196         full_backward_hooks, non_full_backward_hooks = [], []

/home/scemama/QMCTorch/QMCTorch/qmctorch/wavefunction/orbitals/atomic_orbitals.py in forward(self, pos, derivative, sum_grad, sum_hess, one_elec)
    172 
    173         if derivative == [0]:
--> 174             ao = self._compute_ao_values(pos)
    175 
    176         elif derivative == [1]:

/home/scemama/QMCTorch/QMCTorch/qmctorch/wavefunction/orbitals/atomic_orbitals.py in _compute_ao_values(self, pos)
    211 
    212         Y = self.harmonics(xyz)
--> 213         return self._ao_kernel(R, Y)
    214 
    215     def _ao_kernel(self, R, Y):

/home/scemama/QMCTorch/QMCTorch/qmctorch/wavefunction/orbitals/atomic_orbitals.py in _ao_kernel(self, R, Y)
    225         ao = self.norm_cst * R * Y
    226         if self.contract:
--> 227             ao = self._contract(ao)
    228         return ao
    229 

/home/scemama/QMCTorch/QMCTorch/qmctorch/wavefunction/orbitals/atomic_orbitals.py in _contract(self, bas)
    589         """
    590         nbatch = bas.shape[0]
--> 591         bas = self.bas_coeffs * bas
    592         cbas = torch.zeros(nbatch, self.nelec,
    593                            self.norb, device=self.device

RuntimeError: The size of tensor a (49) must match the size of tensor b (41) at non-singleton dimension 2
NicoRenaud commented 8 months ago

@scemama I've finally found time to fix that bug, thanks for flagging it. it should be working now. However i want to refactor the way I read orbital info from pyscf and I also want to support f orbitals that are currently not supported. Thanks !