Added CAS-SCF module

plugins/CASSCF/EZFIO.cfg

@@ -0,0 +1,10 @@
+[energy]
+type: double precision
+doc: "Calculated CAS-SCF energy"
+interface: ezfio
+
+[energy_pt2]
+type: double precision
+doc: "Calculated selected CAS-SCF energy with PT2 correction"
+interface: ezfio
+

plugins/CASSCF/H_apply.irp.f

@@ -0,0 +1,39 @@
+use bitmasks
+BEGIN_SHELL [ /usr/bin/env python ]
+from generate_h_apply import *
+
+s = H_apply("CAS_SD")
+print s
+
+s = H_apply("CAS_SD_selected_no_skip")
+s.set_selection_pt2("epstein_nesbet_2x2")
+s.unset_skip()
+print s
+
+s = H_apply("CAS_SD_selected")
+s.set_selection_pt2("epstein_nesbet_2x2")
+print s
+
+s = H_apply("CAS_SD_PT2")
+s.set_perturbation("epstein_nesbet_2x2")
+print s
+
+
+s = H_apply("CAS_S",do_double_exc=False)
+print s
+
+s = H_apply("CAS_S_selected_no_skip",do_double_exc=False)
+s.set_selection_pt2("epstein_nesbet_2x2")
+s.unset_skip()
+print s
+
+s = H_apply("CAS_S_selected",do_double_exc=False)
+s.set_selection_pt2("epstein_nesbet_2x2")
+print s
+
+s = H_apply("CAS_S_PT2",do_double_exc=False)
+s.set_perturbation("epstein_nesbet_2x2")
+print s
+
+END_SHELL
+

plugins/CASSCF/NEEDED_CHILDREN_MODULES

@@ -0,0 +1 @@
+Generators_CAS Perturbation Selectors_full

plugins/CASSCF/README.rst

@@ -0,0 +1,20 @@
+======
+CASSCF
+======
+
+This module is not a "real" CAS-SCF. It is an orbital optimization step done by :
+
+1) Doing the CAS+SD
+2) Taking one-electron density matrix
+3) Cancelling all active-active rotations
+4) Finding the order which matches with the input MOs
+
+
+Needed Modules
+==============
+.. Do not edit this section It was auto-generated
+.. by the `update_README.py` script.
+Documentation
+=============
+.. Do not edit this section It was auto-generated
+.. by the `update_README.py` script.

plugins/CASSCF/casscf.irp.f

@@ -0,0 +1,220 @@
+program casscf
+  implicit none
+  BEGIN_DOC
+! Optimize MOs and CI coefficients of the CAS
+  END_DOC
+  
+  double precision, allocatable  :: pt2(:), norm_pert(:), H_pert_diag(:)
+  integer(bit_kind), allocatable  :: generators_bitmask_save(:,:,:,:) 
+
+  integer                        :: degree, N_generators_bitmask_save, N_det_ci
+  double precision               :: E_old, E_CI
+  double precision               :: selection_criterion_save, selection_criterion_min_save
+  
+  integer                        :: N_det_old
+  integer                        :: i, j, k, l
+  integer                        :: i_bit, j_bit, i_int, j_int
+  integer(bit_kind), allocatable :: bit_tmp(:), cas_bm(:)
+  character*(64) :: label
+
+  allocate( pt2(N_states), norm_pert(N_states),H_pert_diag(N_states) )
+  allocate( generators_bitmask_save(N_int,2,6,N_generators_bitmask) )
+  allocate( bit_tmp(N_int), cas_bm(N_int) )
+
+  PROVIDE N_det_cas
+  N_det_old = 0
+  pt2 = 1.d0
+  E_CI = 1.d0
+  E_old = 0.d0
+  diag_algorithm = "Lapack"
+  selection_criterion_save = selection_criterion
+  selection_criterion_min_save = selection_criterion_min
+
+
+  cas_bm = 0_bit_kind
+  do i=1,N_cas_bitmask
+    do j=1,N_int
+      cas_bm(j) = ior(cas_bm(j), cas_bitmask(j,1,i))
+      cas_bm(j) = ior(cas_bm(j), cas_bitmask(j,2,i))
+    enddo
+  enddo
+
+  ! Save CAS-SD bitmask
+  generators_bitmask_save = generators_bitmask
+  N_generators_bitmask_save = N_generators_bitmask
+
+  ! Set the CAS bitmask 
+  do i=1,6
+    generators_bitmask(:,:,i,:) = cas_bitmask
+  enddo
+  N_generators_bitmask = N_cas_bitmask
+  SOFT_TOUCH generators_bitmask N_generators_bitmask
+
+
+  ! If the number of dets already in the file is larger than the requested
+  ! number of determinants, truncate the wf
+  if (N_det > N_det_max) then
+    call diagonalize_CI
+    call save_wavefunction
+    psi_det = psi_det_sorted
+    psi_coef = psi_coef_sorted
+    N_det = N_det_max
+    soft_touch N_det psi_det psi_coef
+    call diagonalize_CI
+    call save_wavefunction
+    print *,  'N_det    = ', N_det
+    print *,  'N_states = ', N_states
+    print *,  'PT2      = ', pt2
+    print *,  'E        = ', CI_energy
+    print *,  'E+PT2    = ', CI_energy+pt2
+    print *,  '-----'
+  endif
+
+  ! Start MCSCF iteration
+
+  ! CAS-CI
+  ! ------
+
+  E_old = E_CI
+
+  ! Reset the selection criterion
+  selection_criterion     = selection_criterion_save
+  selection_criterion_min = selection_criterion_min_save
+  SOFT_TOUCH selection_criterion_min selection_criterion selection_criterion_factor
+
+  ! Set the CAS bitmask 
+  do i=1,6
+    generators_bitmask(:,:,i,:) = cas_bitmask
+  enddo
+  N_generators_bitmask = N_cas_bitmask
+  SOFT_TOUCH generators_bitmask N_generators_bitmask
+
+  do while (N_det < N_det_max.and.maxval(abs(pt2(1:N_states))) > pt2_max)
+    N_det_old = N_det
+    call H_apply_CAS_SD_selected_no_skip(pt2, norm_pert, H_pert_diag,  N_states)
+
+    PROVIDE  psi_coef
+    PROVIDE  psi_det
+    PROVIDE  psi_det_sorted
+
+    if (N_det > N_det_max) then
+      psi_det = psi_det_sorted
+      psi_coef = psi_coef_sorted
+      N_det = N_det_max
+      soft_touch N_det psi_det psi_coef
+    endif
+    call diagonalize_CI
+    call save_wavefunction
+    print *,  '======'
+    print *,  'CAS-CI'
+    print *,  '======'
+    print *,  ''
+    print *,  'N_det      = ', N_det
+    print *,  'N_states   = ', N_states
+    print *,  'PT2        = ', pt2
+    print *,  'E(CAS)     = ', CI_energy
+    print *,  'E(CAS)+PT2 = ', CI_energy+pt2
+    print *,  '-----'
+    print *,  ''
+    E_CI = sum(CI_energy(1:N_states)+pt2(1:N_states))/dble(N_states)
+
+    call ezfio_set_casscf_energy(CI_energy(1))
+    if (abort_all) then
+      exit
+    endif
+    if (N_det == N_det_old) then
+      exit
+    endif
+
+  enddo
+
+  ! Super-CI
+  ! --------
+
+  selection_criterion_min = 1.d-12
+  selection_criterion     = 1.d-12
+
+  ! Set the CAS bitmask
+  generators_bitmask = generators_bitmask_save
+  N_generators_bitmask = N_generators_bitmask_save
+  SOFT_TOUCH generators_bitmask N_generators_bitmask selection_criterion selection_criterion_min selection_criterion_factor
+
+  N_det_ci = N_det
+
+  call H_apply_CAS_SD_selected(pt2, norm_pert, H_pert_diag,  N_states)
+  
+
+  do i=1,mo_tot_num
+   i_int = ishft(i-1,-bit_kind_shift)+1
+   i_bit = j-ishft(i_int-1,bit_kind_shift)-1
+   bit_tmp(:) = 0_bit_kind
+   bit_tmp(i_int) = ibset(0_bit_kind,i_bit)
+   if (iand(bit_tmp(i_int), cas_bm(i_int)) == 0_bit_kind) then
+     ! Not a CAS MO
+     cycle
+   endif
+
+   do j=1,mo_tot_num
+     if (j == i) then
+       cycle
+     endif
+     j_int = ishft(j-1,-bit_kind_shift)+1
+     j_bit = j-ishft(j_int-1,bit_kind_shift)-1
+     bit_tmp(:) = 0_bit_kind
+     bit_tmp(j_int) = ibset(0_bit_kind,j_bit)
+     if (iand(bit_tmp(j_int), cas_bm(j_int)) == 0_bit_kind) then
+       ! Not a CAS MO
+       cycle
+     endif
+     ! Now, both i and j are MOs of the CAS. De-couple them in the DM
+     one_body_dm_mo(i,j) = 0.d0
+   enddo
+
+  enddo
+
+  SOFT_TOUCH one_body_dm_mo
+
+  double precision               :: mx, ov
+  double precision, allocatable  :: mo_coef_old(:,:)
+  integer, allocatable           :: iorder(:)
+  logical, allocatable           :: selected(:)
+  allocate( mo_coef_old(size(mo_coef,1), size(mo_coef,2)), iorder(mo_tot_num), selected(mo_tot_num) )
+  mo_coef_old = mo_coef
+  label = "Canonical"
+  call mo_as_eigvectors_of_mo_matrix(one_body_dm_mo,size(one_body_dm_mo,1),size(one_body_dm_mo,2),label,-1)
+  selected = .False.
+  do j=1,mo_tot_num
+   mx = -1.d0
+   iorder(j) = j
+   do i=1,mo_tot_num
+    if (selected(i)) then
+      cycle
+    endif
+    ov = 0.d0
+    do l=1,ao_num
+     do k=1,ao_num
+      ov = ov + mo_coef_old(k,j) * ao_overlap(k,l) * mo_coef(l,i)
+     enddo
+    enddo
+    ov= dabs(ov)
+    if (ov > mx) then
+      mx = ov
+      iorder(j) = i
+    endif
+   enddo
+   selected( iorder(j) ) = .True.
+  enddo
+  mo_coef_old = mo_coef
+  do i=1,mo_tot_num
+    mo_coef(:,i) = mo_coef_old(:,iorder(i))
+  enddo
+
+  call save_mos
+
+  call write_double(6,E_CI,"Energy(CAS)")
+
+  deallocate( mo_coef_old )
+  deallocate( pt2, norm_pert,H_pert_diag )
+  deallocate( generators_bitmask_save )
+  deallocate( bit_tmp, cas_bm, iorder )
+end

src/Bitmask/bitmasks.irp.f

@@ -262,13 +262,7 @@ END_PROVIDER
  logical                        :: exists
  integer                        :: j,i
  integer :: i_hole,i_part,i_gen
- PROVIDE ezfio_filename
+
-!do j = 1, N_int
-! inact_bitmask(j,1) = xor(generators_bitmask(j,1,1,1),cas_bitmask(j,1,1))
-! inact_bitmask(j,2) = xor(generators_bitmask(j,2,1,1),cas_bitmask(j,2,1))
-! virt_bitmask(j,1) = xor(generators_bitmask(j,1,2,1),cas_bitmask(j,1,1))
-! virt_bitmask(j,2) = xor(generators_bitmask(j,2,2,1),cas_bitmask(j,2,1))
-!enddo
  n_inact_orb = 0
  n_virt_orb = 0
  if(N_generators_bitmask == 1)then

src/Integrals_Bielec/map_integrals.irp.f

@@ -424,15 +424,6 @@ integer*8 function get_mo_map_size()
   get_mo_map_size = mo_integrals_map % n_elements
 end
 
-subroutine clear_mo_map
-  implicit none
-  BEGIN_DOC
-  ! Frees the memory of the MO map
-  END_DOC
-  call map_deinit(mo_integrals_map)
-  FREE mo_integrals_map
-end
-
 BEGIN_TEMPLATE
 
 subroutine dump_$ao_integrals(filename)

src/Integrals_Bielec/mo_bi_integrals.irp.f

@@ -503,3 +503,14 @@ BEGIN_PROVIDER [ double precision, mo_bielec_integral_schwartz,(mo_tot_num,mo_to
   enddo
   
 END_PROVIDER
+
+
+subroutine clear_mo_map
+  implicit none
+  BEGIN_DOC
+  ! Frees the memory of the MO map
+  END_DOC
+  call map_deinit(mo_integrals_map)
+  FREE mo_integrals_map mo_bielec_integral_schwartz mo_bielec_integral_jj mo_bielec_integral_jj_anti
+  FREE mo_bielec_integral_jj_exchange mo_bielec_integrals_in_map
+end

src/MO_Basis/utils.irp.f

@@ -98,7 +98,6 @@ subroutine mo_as_eigvectors_of_mo_matrix(matrix,n,m,label,sign)
   call write_time(output_mo_basis)
   
   mo_label = label
-  SOFT_TOUCH mo_coef mo_label
 end
 
 subroutine mo_as_eigvectors_of_mo_matrix_sort_by_observable(matrix,observable,n,m,label)