MCSCF density

EZFIO.tar.gz

@@ -0,0 +1 @@
+EZFIO.1.0.7.tar.gz

bin/to_ezfio.py

@@ -137,6 +137,7 @@ def write_ezfioFile(res,filename):
   ezfio.mo_basis_mo_coef = MoMatrix
 
 # Determinants
+  det_thr = 1.e-16
   closed_mos = res.closed_mos
   nactive = ezfio.get_mo_basis_mo_active_num()
   dets_a = []
@@ -150,7 +151,9 @@ def write_ezfioFile(res,filename):
     for x in d['beta']:
       if x not in closed_mos:
         dnew_b.append(x+1)
-    for x in range(nactive-len(dnew_b)):
+    for x in range(res.num_alpha-len(dnew_a)-len(closed_mos)):
+      dnew_a.append(0)
+    for x in range(res.num_alpha-len(dnew_b)-len(closed_mos)):
       dnew_b.append(0)
     dets_a.append( dnew_a )
     dets_b.append( dnew_b )
@@ -158,6 +161,12 @@ def write_ezfioFile(res,filename):
   coef = reduce(lambda x, y: x+y,res.det_coefficients,[])
 
   if len(dets_a[0]) > 0:
+   for i in xrange(len(coef),0,-1):
+     i -= 1
+     if coef[i] == 0.:
+       dets_a.pop(i)
+       dets_b.pop(i)
+       coef.pop(i)
    ezfio.determinants_det_num = len(dets_a)
    ezfio.determinants_det_coef = coef
    ezfio.determinants_det_occ = dets_a+dets_b

eplf.config

@@ -10,7 +10,7 @@ nuclei
 
 determinants
   det_num                integer
-  det_occ                integer           (mo_basis_mo_active_num,determinants_det_num,2)
+  det_occ                integer           (electrons_elec_alpha_num-mo_basis_mo_closed_num,determinants_det_num,2)
   det_coef               real              (determinants_det_num)
 
 ao_basis

src/density.irp.f

@@ -16,11 +16,9 @@ BEGIN_PROVIDER [ double precision, density_grad_p, (3) ]
 ! Gradient of the density at the current point 
  END_DOC
 
- integer :: i, l
+ density_grad_p(1) = density_alpha_grad_p(1) + density_beta_grad_p(1)
-
+ density_grad_p(2) = density_alpha_grad_p(2) + density_beta_grad_p(2)
- do l=1,3
+ density_grad_p(3) = density_alpha_grad_p(3) + density_beta_grad_p(3)
-  density_grad_p(l) = density_alpha_grad_p(l) + density_beta_grad_p(l)
- enddo
 
 END_PROVIDER
 
@@ -31,10 +29,25 @@ BEGIN_PROVIDER [ real, density_alpha_value_p ]
 ! Value of the alpha density at the current point 
  END_DOC
 
- density_alpha_value_p = 0.
  integer :: i
- do i=1,elec_alpha_num
+
-  density_alpha_value_p = density_alpha_value_p + mo_value_p(i)**2
+ density_alpha_value_p = 0.
+
+ do i=1,mo_closed_num
+  density_alpha_value_p += mo_value_p(i)**2
+ enddo
+
+! TODO vectorization
+ integer :: k,j,l
+ real    :: buffer
+ do k=1,det_num
+  do l=1,det_num
+   buffer = 0.
+   do i=1,mo_active_num
+    buffer += mo_value_p(det(i,k,1))*mo_value_p(det(i,l,1))
+   enddo
+   density_alpha_value_p += det_coef(k)*det_coef(l)*buffer
+  enddo
  enddo
 
 END_PROVIDER
@@ -47,8 +60,21 @@ BEGIN_PROVIDER [ real, density_beta_value_p ]
 
  density_beta_value_p = 0.
  integer :: i
- do i=1,elec_beta_num
+ do i=1,mo_closed_num
-  density_beta_value_p = density_beta_value_p + mo_value_p(i)**2
+  density_beta_value_p += mo_value_p(i)**2
+ enddo
+
+! TODO vectorization
+ integer :: k,j,l
+ real    :: buffer
+ do k=1,det_num
+  do l=1,det_num
+   buffer = 0.
+   do i=1,mo_active_num
+    buffer += mo_value_p(det(i,k,2))*mo_value_p(det(i,l,2))
+   enddo
+   density_beta_value_p += det_coef(k)*det_coef(l)*buffer
+  enddo
  enddo
 
 END_PROVIDER
@@ -60,15 +86,37 @@ BEGIN_PROVIDER [ double precision, density_alpha_grad_p, (3) ]
 ! Gradient of the density at the current point 
  END_DOC
 
- integer :: i, l
+ integer :: i
 
- do l=1,3
+ density_alpha_grad_p(1) = 0.
-  density_alpha_grad_p(l) = 0.
+ density_alpha_grad_p(2) = 0.
+ density_alpha_grad_p(3) = 0.
+
+ do i=1,mo_closed_num
+   density_alpha_grad_p(1) += 2.*mo_grad_p(i,1)*mo_value_p(i)
+   density_alpha_grad_p(2) += 2.*mo_grad_p(i,2)*mo_value_p(i)
+   density_alpha_grad_p(3) += 2.*mo_grad_p(i,3)*mo_value_p(i)
  enddo
 
- do i=1,elec_alpha_num
+! TODO vectorization
-  do l=1,3
+ integer :: k,j,l
-   density_alpha_grad_p(l) = density_alpha_grad_p(l) + 2.*mo_grad_p(i,l)*mo_value_p(i)
+ real    :: buffer(3)
+ do k=1,det_num
+  do l=1,det_num
+   buffer(1) = 0.
+   buffer(2) = 0.
+   buffer(3) = 0.
+   do i=1,mo_active_num
+    buffer(1) += mo_grad_p(det(i,k,1),1)*mo_value_p(det(i,l,1)) &
+               + mo_grad_p(det(i,l,1),1)*mo_value_p(det(i,k,1)) 
+    buffer(2) += mo_grad_p(det(i,k,1),2)*mo_value_p(det(i,l,1)) &
+               + mo_grad_p(det(i,l,1),2)*mo_value_p(det(i,k,1)) 
+    buffer(3) += mo_grad_p(det(i,k,1),3)*mo_value_p(det(i,l,1)) &
+               + mo_grad_p(det(i,l,1),3)*mo_value_p(det(i,k,1)) 
+   enddo
+   density_alpha_grad_p(1) += det_coef(k)*det_coef(l)*buffer(1)
+   density_alpha_grad_p(2) += det_coef(k)*det_coef(l)*buffer(2)
+   density_alpha_grad_p(3) += det_coef(k)*det_coef(l)*buffer(3)
   enddo
  enddo
 
@@ -81,15 +129,37 @@ BEGIN_PROVIDER [ double precision, density_beta_grad_p, (3) ]
 ! Gradient of the density at the current point 
  END_DOC
 
- integer :: i, l
+ integer :: i
 
- do l=1,3
+ density_beta_grad_p(1) = 0.
-  density_beta_grad_p(l) = 0.
+ density_beta_grad_p(2) = 0.
+ density_beta_grad_p(3) = 0.
+
+ do i=1,mo_closed_num
+   density_beta_grad_p(1) += 2.*mo_grad_p(i,1)*mo_value_p(i)
+   density_beta_grad_p(2) += 2.*mo_grad_p(i,2)*mo_value_p(i)
+   density_beta_grad_p(3) += 2.*mo_grad_p(i,3)*mo_value_p(i)
  enddo
 
- do i=1,elec_beta_num
+! TODO vectorization
-  do l=1,3
+ integer :: k,j,l
-   density_beta_grad_p(l) = density_beta_grad_p(l) + 2.*mo_grad_p(i,l)*mo_value_p(i)
+ real    :: buffer(3)
+ do k=1,det_num
+  do l=1,det_num
+   buffer(1) = 0.
+   buffer(2) = 0.
+   buffer(3) = 0.
+   do i=1,mo_active_num
+    buffer(1) += mo_grad_p(det(i,k,2),1)*mo_value_p(det(i,l,2)) &
+               + mo_grad_p(det(i,l,2),1)*mo_value_p(det(i,k,2)) 
+    buffer(2) += mo_grad_p(det(i,k,2),2)*mo_value_p(det(i,l,2)) &
+               + mo_grad_p(det(i,l,2),2)*mo_value_p(det(i,k,2)) 
+    buffer(3) += mo_grad_p(det(i,k,2),3)*mo_value_p(det(i,l,2)) &
+               + mo_grad_p(det(i,l,2),3)*mo_value_p(det(i,k,2)) 
+   enddo
+   density_beta_grad_p(1) += det_coef(k)*det_coef(l)*buffer(1)
+   density_beta_grad_p(2) += det_coef(k)*det_coef(l)*buffer(2)
+   density_beta_grad_p(3) += det_coef(k)*det_coef(l)*buffer(3)
   enddo
  enddo
 
@@ -104,17 +174,34 @@ BEGIN_PROVIDER [ double precision, density_alpha_lapl_p ]
 ! Laplacian of the density at the current point 
  END_DOC
 
- integer :: i, l
+ integer :: i
 
  density_alpha_lapl_p = 0.
 
- do i=1,elec_alpha_num
+ do i=1,mo_closed_num
-  do l=1,3
+  density_alpha_lapl_p += mo_grad_p(i,1)**2
-   density_alpha_lapl_p = density_alpha_lapl_p + mo_grad_p(i,l)**2
+  density_alpha_lapl_p += mo_grad_p(i,2)**2
-  enddo
+  density_alpha_lapl_p += mo_grad_p(i,3)**2
-  density_alpha_lapl_p = density_alpha_lapl_p + mo_value_p(i)*mo_lapl_p(i)
+  density_alpha_lapl_p += mo_value_p(i)*mo_lapl_p(i)
+ enddo
+ density_alpha_lapl_p *= 2.
+
+! TODO vectorization
+ integer :: k,j,l
+ real    :: buffer
+ do k=1,det_num
+  do l=1,det_num
+   buffer = 0.
+   do i=1,mo_active_num
+    buffer += 2.*(mo_grad_p(det(i,k,1),1)*mo_grad_p(det(i,l,1),1) &
+                + mo_grad_p(det(i,k,1),2)*mo_grad_p(det(i,l,1),2) &
+                + mo_grad_p(det(i,k,1),3)*mo_grad_p(det(i,l,1),3))&
+           + mo_value_p(det(i,k,1))*mo_lapl_p(det(i,l,1)) &
+           + mo_value_p(det(i,l,1))*mo_lapl_p(det(i,k,1)) 
+   enddo
+   density_alpha_lapl_p += det_coef(k)*det_coef(l)*buffer
+  enddo
  enddo
- density_alpha_lapl_p = 2.*density_alpha_lapl_p
 
 END_PROVIDER
 
@@ -125,17 +212,34 @@ BEGIN_PROVIDER [ double precision, density_beta_lapl_p ]
 ! Laplacian of the density at the current point 
  END_DOC
 
- integer :: i, l
+ integer :: i
 
  density_beta_lapl_p = 0.
 
- do i=1,elec_beta_num
+ do i=1,mo_closed_num
-  do l=1,3
+  density_beta_lapl_p += mo_grad_p(i,1)**2
-   density_beta_lapl_p = density_beta_lapl_p + mo_grad_p(i,l)**2
+  density_beta_lapl_p += mo_grad_p(i,2)**2
-  enddo
+  density_beta_lapl_p += mo_grad_p(i,3)**2
-  density_beta_lapl_p = density_beta_lapl_p + mo_value_p(i)*mo_lapl_p(i)
+  density_beta_lapl_p += mo_value_p(i)*mo_lapl_p(i)
+ enddo
+ density_beta_lapl_p *= 2.
+
+! TODO vectorization
+ integer :: k,j,l
+ real    :: buffer
+ do k=1,det_num
+  do l=1,det_num
+   buffer = 0.
+   do i=1,mo_active_num
+    buffer += 2.*(mo_grad_p(det(i,k,2),1)*mo_grad_p(det(i,l,2),1) &
+                + mo_grad_p(det(i,k,2),2)*mo_grad_p(det(i,l,2),2) &
+                + mo_grad_p(det(i,k,2),3)*mo_grad_p(det(i,l,2),3))&
+           + mo_value_p(det(i,k,2))*mo_lapl_p(det(i,l,2)) &
+           + mo_value_p(det(i,l,2))*mo_lapl_p(det(i,k,2)) 
+   enddo
+   density_beta_lapl_p += det_coef(k)*det_coef(l)*buffer
+  enddo
  enddo
- density_beta_lapl_p = 2.*density_beta_lapl_p
 
 END_PROVIDER
 
@@ -145,8 +249,6 @@ BEGIN_PROVIDER [ double precision, density_lapl_p ]
 ! Laplacian of the density at the current point 
  END_DOC
 
- integer :: i, l
-
  density_lapl_p = density_alpha_lapl_p + density_beta_lapl_p
 
 END_PROVIDER

src/det.irp.f

@@ -18,12 +18,16 @@ BEGIN_PROVIDER  [ integer, det, (elec_alpha_num-mo_closed_num,det_num,2) ]
 
   BEGIN_DOC
 ! det : Description of the active orbitals of the determinants
+
 ! det_coef  : Determinant coefficients
   END_DOC
 
- det = 0
+ if (elec_alpha_num > mo_closed_num) then
+  det = 0
+  call get_determinants_det_occ(det)
+ endif
  det_coef = 0.
- call get_determinants_det_occ(det)
+ det_coef(1) = 1.
  call get_determinants_det_coef(det_coef)
 
 END_PROVIDER