Created
August 24, 2011 08:08
-
-
Save gabrielelanaro/1167545 to your computer and use it in GitHub Desktop.
mp2 impl
This file contains hidden or bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
| from calculator_hf import HartreeFock | |
| import numpy as np | |
| class MP2: | |
| def __init__(self, scf, text, parameters): | |
| """ Takes as input *scf* that is a converged SCF and the usual | |
| output function *text* and the parameter dictionary | |
| *parameters. | |
| """ | |
| self.text = text | |
| self.parameters = parameters | |
| self.scf = scf | |
| def get_correction(self): | |
| # Number of orbitals | |
| occ = self.scf.calculator.occupied | |
| virt = self.scf.calculator.getUnoccupied() | |
| tot = occ+virt | |
| en = self.scf.orbitals.energiesAlpha | |
| res = 0 | |
| epairs = np.zeros((occ, occ), "d") | |
| for i in range(occ): | |
| for j in range(occ): | |
| for a in range(occ, tot): | |
| for b in range(occ, tot): | |
| iajb = self.ijkl(i,a,j,b) | |
| ibja = self.ijkl(i,b,j,a) | |
| epairs[i,j] += iajb*(2*iajb - ibja) / (en[i] + en[j] -en[a] -en[b]) | |
| return sum(sum(epairs)) | |
| def four_index_transformation(self): | |
| """ I have to form the tensor (ij|kl). | |
| It is in chemists' notation | |
| """ | |
| scf = self.scf | |
| # Orbital coefficients | |
| orbs = scf.orbitals.coefficientsAlpha.transpose() | |
| eri = scf.calculator.two.V | |
| nbas = orbs.shape[0] | |
| def iqrs(i, q, r, s): | |
| return sum(C_p*eri[p,q,r,s] for p, C_p in enumerate(orbs[i])) | |
| def ijrs(i, j, r, s): | |
| return sum(C_q*iqrs(i,q,r,s) for q, C_q in enumerate(orbs[j])) | |
| def ijks(i,j,k,s): | |
| return sum(C_r*ijrs(i,j,r,s) for r, C_r in enumerate(orbs[k])) | |
| def ijkl(i,j,k,l): | |
| return sum(C_s*ijks(i,j,k,s) for s, C_s in enumerate(orbs[l])) | |
| #self.ijkl = np.fromfunction(ijkl, (nbas, nbas, nbas, nbas), dtype=int) | |
| self.ijkl = ijkl |
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment