Permalink
Browse files

Implement Gaussian::Primitive#overlap

  • Loading branch information...
1 parent 57cc455 commit 97aa0663c2ca0fe4e58a47e2aebda87a3e3e8124 @ktns committed May 12, 2012
Showing with 38 additions and 2 deletions.
  1. +38 −2 lib/ruphy/orbital/gaussian.rb
@@ -29,8 +29,44 @@ def normalization_factor
overlap(self)**(-0.5)
end
- def overlap other
- raise NotImplementedError
+ def overlap o
+ z=@zeta+o.zeta
+ c=(@center*@zeta+o.center*o.zeta)*z**-1.0
+ [@center-c,o.center-c,@momenta,o.momenta].map(&:to_a).transpose.map do |a,b,m,n|
+ integral(a,b,m,n,z)
+ end.reduce(:*) * exp(@zeta*o.zeta/z*(@center-o.center).r**2)
+ end
+
+ include Math
+ private
+ # calculate one dimensional integral over entire space of
+ # (x-a)^m(x-b)^n exp(-z*x^2)
+ def integral a,b,m,n,z
+ case m
+ when 0
+ case n
+ when 0
+ sqrt(PI/z)
+ when 1
+ -b*sqrt(PI/z)
+ else
+ (b**2+(n-1)/2.0/z) * integral(a,b,m,n-2,z) -
+ 2*b * integral(a,b,m,n-1,z)
+ end
+ when 1
+ case n
+ when 0
+ -a*sqrt(PI/z)
+ when 1
+ (a*b+0.5/z)*sqrt(PI/z)
+ else
+ (b**2+(n-1)/2.0/z) * integral(a,b,m,n-2,z) -
+ 2*b * integral(a,b,m,n-1,z)
+ end
+ else
+ (a**2+(m-1)/2.0/z) * integral(a,b,m-2,n,z) -
+ 2*a * integral(a,b,m-1,n,z)
+ end
end
end

0 comments on commit 97aa066

Please sign in to comment.