Merge branch 'dev-stable' of github.com:QuantumPackage/qp2 into dev-stable

src/casscf_cipsi/README.rst

@@ -4,13 +4,15 @@ casscf
 
 |CASSCF| program with the CIPSI algorithm.
 
-Example of inputs
------------------
+
+Example of inputs for GROUND STATE calculations
+-----------------------------------------------
+NOTICE :: FOR EXCITED STATES CALCULATIONS SEE THE FILE "example_casscf_multistate.sh"
 
 a) Small active space : standard CASSCF 
 ---------------------------------------
 Let's do O2 (triplet) in aug-cc-pvdz with the following geometry (xyz format, Bohr units)
-3
+2
 
  O           0.0000000000        0.0000000000       -1.1408000000
  O           0.0000000000        0.0000000000        1.1408000000
@@ -45,3 +47,4 @@ qp set casscf_cipsi small_active_space False
 qp run casscf | tee ${EZFIO_FILE}.casscf_large.out
 # you should find around -149.9046
 
+

src/casscf_cipsi/casscf.irp.f

@@ -54,14 +54,24 @@ subroutine run
 
     call write_time(6)
     call write_int(6,iteration,'CAS-SCF iteration = ')
-    call write_double(6,energy,'CAS-SCF energy = ')
+    call write_double(6,energy,'State-average CAS-SCF energy = ')
 !    if(n_states == 1)then
 !     call ezfio_get_casscf_cipsi_energy_pt2(E_PT2)
 !     call ezfio_get_casscf_cipsi_energy(PT2)
+     double precision :: delta_E_istate, e_av
+     e_av = 0.d0
      do istate=1,N_states
-     call write_double(6,E_PT2(istate),'E + PT2 energy = ')
-     call write_double(6,PT2(istate),'  PT2          = ')
+      e_av += state_average_weight(istate) * Ev(istate) 
+      if(istate.gt.1)then
+       delta_E_istate = E_PT2(istate)  - E_PT2(1)
+      write(*,'(A6,I2,A18,F16.10)')'state ',istate,' Delta E+PT2    = ',delta_E_istate
+      endif
+      write(*,'(A6,I2,A18,F16.10)')'state ',istate,' E + PT2 energy = ',E_PT2(istate)
+      write(*,'(A6,I2,A18,F16.10)')'state ',istate,'     PT2 energy = ',PT2(istate)
+!      call write_double(6,E_PT2(istate),'E + PT2 energy = ')
+!      call write_double(6,PT2(istate),'  PT2          = ')
      enddo
+     call write_double(6,e_av,'State-average CAS-SCF energy bis = ')
      call write_double(6,pt2_max,' PT2_MAX       = ')
 !    endif
 
@@ -99,8 +109,8 @@ subroutine run
 
     mo_coef = NewOrbs
     mo_occ  = occnum
-    call save_mos
     if(.not.converged)then
+     call save_mos
      iteration += 1
      if(norm_grad_vec2.gt.0.01d0)then
       N_det = N_states

src/casscf_cipsi/example_casscf_multistate.sh

@@ -0,0 +1,66 @@
+# This is an example for MULTI STATE CALCULATION STATE AVERAGE CASSCF 
+# We will compute 3 states on the O2 molecule 
+# The Ground state and 2 degenerate excited states 
+# Please follow carefully the tuto :) 
+
+##### PREPARING THE EZFIO 
+# Set the path to your QP2 directory 
+QP_ROOT=my_fancy_path 
+source ${QP_ROOT}/quantum_package.rc 
+# Create the EZFIO folder
+qp create_ezfio -b aug-cc-pvdz O2.xyz -m 3 -a -o O2_avdz_multi_state
+# Start with ROHF orbitals 
+qp run scf # ROHF energy : -149.619992871398
+# Freeze the 1s orbitals of the two oxygen
+qp set_frozen_core 
+
+##### PREPARING THE ORBITALS WITH NATURAL ORBITALS OF A CIS 
+# Tell that you want 3 states in your WF
+qp set determinants n_states 3 
+# Run a CIS wave function to start your calculation 
+qp run cis | tee ${EZFIO_FILE}.cis_3_states.out # -149.6652601409258 -149.4714726176746 -149.4686165431939
+# Save the STATE AVERAGE natural orbitals for having a balanced description 
+# This will also order the orbitals according to their occupation number 
+# Which makes the active space selection easyer !
+qp run save_natorb | tee ${EZFIO_FILE}.natorb_3states.out 
+
+##### PREPARING A CIS GUESS WITHIN THE ACTIVE SPACE 
+# Set an active space which has the most of important excitations 
+# and that maintains symmetry : the ACTIVE ORBITALS are from """6 to 13"""
+
+# YOU FIRST FREEZE THE VIRTUALS THAT ARE NOT IN THE ACTIVE SPACE 
+# !!!!! WE SET TO "-D" for DELETED !!!!
+qp set_mo_class -c "[1-5]" -a "[6-13]" -d "[14-46]" 
+# You create a guess of CIS type WITHIN THE ACTIVE SPACE 
+qp run cis | tee ${EZFIO_FILE}.cis_3_states_active_space.out # -149.6515472533511 -149.4622878024821 -149.4622878024817
+# You tell to read the WFT stored (i.e. the guess we just created) 
+qp set determinants read_wf True
+
+##### DOING THE CASSCF 
+### SETTING PROPERLY THE ACTIVE SPACE FOR CASSCF 
+# You set the active space WITH THE VIRTUAL ORBITALS !!!
+# !!!!! NOW WE SET TO "-v" for VIRTUALS  !!!!!
+qp set_mo_class -c "[1-5]" -a "[6-13]" -v "[14-46]" 
+
+# You tell that it is a small actice space so the CIPSI can take all Slater determinants 
+qp set casscf_cipsi small_active_space True 
+# You specify the output file 
+output=${EZFIO_FILE}.casscf_3states.out 
+# You run the CASSCF calculation 
+qp run casscf  | tee ${output} # -149.7175867510 -149.5059010227 -149.5059010226
+
+# Some grep in order to get some numbers useful to check convergence 
+# State average energy 
+grep "State-average CAS-SCF energy =" $output | cut -d "=" -f 2 > data_e_average
+# Delta E anticipated for State-average energy, only usefull to check convergence
+grep "Predicted energy improvement =" $output | cut -d "=" -f 2 > data_improve
+# Ground state energy
+grep "state  1 E + PT2 energy"        $output | cut -d "=" -f 2 > data_1
+# First excited state energy
+grep "state  2 E + PT2 energy"        $output | cut -d "=" -f 2 > data_2
+# First excitation energy 
+grep "state  2 Delta E+PT2"           $output | cut -d "=" -f 2 > data_delta_E2
+# Second excited state energy
+grep "state  3 E + PT2 energy"        $output | cut -d "=" -f 2  > data_3
+# Second excitation energy 
+grep "state  3 Delta E+PT2"           $output | cut -d "=" -f 2 > data_delta_E3

src/casscf_cipsi/neworbs.irp.f

@@ -226,27 +226,28 @@ BEGIN_PROVIDER [real*8, Umat, (mo_num,mo_num) ]
   end do
   
   ! Form the exponential
+  call exp_matrix_taylor(Tmat,mo_num,Umat,converged)
 
-  Tpotmat(:,:)=0.D0
-  Umat(:,:)   =0.D0
-  do i=1,mo_num
-    Tpotmat(i,i)=1.D0
-    Umat(i,i)   =1.d0
-  end do
-  iter=0
-  converged=.false.
-  do while (.not.converged)
-    iter+=1
-    f = 1.d0 / dble(iter)
-    Tpotmat2(:,:) = Tpotmat(:,:) * f
-    call dgemm('N','N', mo_num,mo_num,mo_num,1.d0,                   &
-        Tpotmat2, size(Tpotmat2,1),                                  &
-        Tmat, size(Tmat,1), 0.d0,                                    &
-        Tpotmat, size(Tpotmat,1))
-    Umat(:,:) = Umat(:,:) + Tpotmat(:,:)
-    
-    converged = ( sum(abs(Tpotmat(:,:))) < 1.d-6).or.(iter>30)
-  end do
+!  Tpotmat(:,:)=0.D0
+!  Umat(:,:)   =0.D0
+!  do i=1,mo_num
+!    Tpotmat(i,i)=1.D0
+!    Umat(i,i)   =1.d0
+!  end do
+!  iter=0
+!  converged=.false.
+!  do while (.not.converged)
+!    iter+=1
+!    f = 1.d0 / dble(iter)
+!    Tpotmat2(:,:) = Tpotmat(:,:) * f
+!    call dgemm('N','N', mo_num,mo_num,mo_num,1.d0,                   &
+!        Tpotmat2, size(Tpotmat2,1),                                  &
+!        Tmat, size(Tmat,1), 0.d0,                                    &
+!        Tpotmat, size(Tpotmat,1))
+!    Umat(:,:) = Umat(:,:) + Tpotmat(:,:)
+!    
+!    converged = ( sum(abs(Tpotmat(:,:))) < 1.d-6).or.(iter>30)
+!  end do
 END_PROVIDER
   
 

src/two_body_rdm/act_2_transition_rdm.irp.f

@@ -1,9 +1,9 @@
  BEGIN_PROVIDER [double precision, act_2_rdm_trans_spin_trace_mo, (n_act_orb,n_act_orb,n_act_orb,n_act_orb,N_states,N_states)]
  implicit none
  BEGIN_DOC
-! act_2_rdm_trans_spin_trace_mo(i,j,k,l,istate) =  STATE SPECIFIC physicist notation for 2rdm_trans 
+! act_2_rdm_trans_spin_trace_mo(i,j,k,l,istate,jstate) =  STATE SPECIFIC physicist notation for 2rdm_trans 
 ! 
-! \sum_{\sigma,\sigma'}<Psi_{istate}| a^{\dagger}_{i \sigma} a^{\dagger}_{j \sigma'} a_{l \sigma'} a_{k \sigma} |Psi_{istate}>
+! \sum_{\sigma,\sigma'}<Psi_{istate}| a^{\dagger}_{i \sigma} a^{\dagger}_{j \sigma'} a_{l \sigma'} a_{k \sigma} |Psi_{jstate}>
 !
 ! WHERE ALL ORBITALS (i,j,k,l) BELONGS TO AN ACTIVE SPACE DEFINED BY "list_act" 
 !

src/utils/linear_algebra.irp.f

@@ -1897,3 +1897,140 @@ end do
 
 end subroutine pivoted_cholesky
 
+subroutine exp_matrix(X,n,exp_X)
+ implicit none
+ double precision, intent(in) :: X(n,n)
+ integer, intent(in):: n
+ double precision, intent(out):: exp_X(n,n)
+ BEGIN_DOC
+ ! exponential of the matrix X: X has to be ANTI HERMITIAN !! 
+ ! 
+ ! taken from Hellgaker, jorgensen, Olsen book
+ !
+ ! section evaluation of matrix exponential (Eqs. 3.1.29 to 3.1.31)
+ END_DOC
+ integer :: i
+ double precision, allocatable :: r2_mat(:,:),eigvalues(:),eigvectors(:,:)
+ double precision, allocatable :: matrix_tmp1(:,:),eigvalues_mat(:,:),matrix_tmp2(:,:)
+ include 'constants.include.F'
+ allocate(r2_mat(n,n),eigvalues(n),eigvectors(n,n))
+ allocate(eigvalues_mat(n,n),matrix_tmp1(n,n),matrix_tmp2(n,n))
+
+ ! r2_mat = X^2 in the 3.1.30 
+ call get_A_squared(X,n,r2_mat)
+ call lapack_diagd(eigvalues,eigvectors,r2_mat,n,n) 
+ eigvalues=-eigvalues
+
+ if(.False.)then
+ !!! For debugging and following the book intermediate
+   ! rebuilding the matrix : X^2 = -W t^2 W^T as in 3.1.30 
+   ! matrix_tmp1 = W t^2 
+   print*,'eigvalues = '
+   do i = 1, n
+    print*,i,eigvalues(i)
+    write(*,'(100(F16.10,X))')eigvectors(:,i)
+   enddo
+   eigvalues_mat=0.d0
+   do i = 1,n
+    ! t = dsqrt(t^2) where t^2 are eigenvalues of X^2
+    eigvalues(i) = dsqrt(eigvalues(i)) 
+    eigvalues_mat(i,i) = eigvalues(i)*eigvalues(i)
+   enddo
+   call dgemm('N','N',n,n,n,1.d0,eigvectors,size(eigvectors,1), &
+   eigvalues_mat,size(eigvalues_mat,1),0.d0,matrix_tmp1,size(matrix_tmp1,1))
+   call dgemm('N','T',n,n,n,-1.d0,matrix_tmp1,size(matrix_tmp1,1), &
+   eigvectors,size(eigvectors,1),0.d0,matrix_tmp2,size(matrix_tmp2,1))
+   print*,'r2_mat new = '
+   do i = 1, n
+    write(*,'(100(F16.10,X))')matrix_tmp2(:,i)
+   enddo
+ endif
+
+ ! building the exponential 
+ ! exp(X) = W cos(t) W^T + W t^-1 sin(t) W^T X as in Eq. 3.1.31
+ ! matrix_tmp1 = W cos(t) 
+ do i = 1,n
+  eigvalues_mat(i,i) = dcos(eigvalues(i))
+ enddo
+ ! matrix_tmp2 = W cos(t)
+ call dgemm('N','N',n,n,n,1.d0,eigvectors,size(eigvectors,1), &
+ eigvalues_mat,size(eigvalues_mat,1),0.d0,matrix_tmp1,size(matrix_tmp1,1))
+ ! matrix_tmp2 = W cos(t) W^T
+ call dgemm('N','T',n,n,n,-1.d0,matrix_tmp1,size(matrix_tmp1,1), &
+ eigvectors,size(eigvectors,1),0.d0,matrix_tmp2,size(matrix_tmp2,1))
+ exp_X = matrix_tmp2
+ ! matrix_tmp2 = W t^-1 sin(t) W^T X
+ do i = 1,n
+  if(dabs(eigvalues(i)).gt.1.d-4)then
+   eigvalues_mat(i,i) = dsin(eigvalues(i))/eigvalues(i)
+  else ! Taylor development of sin(x)/x near x=0 = 1 - x^2/6
+   eigvalues_mat(i,i) = 1.d0 - eigvalues(i)*eigvalues(i)*c_1_3*0.5d0
+  endif
+ enddo
+ ! matrix_tmp1 = W t^-1 sin(t) 
+ call dgemm('N','N',n,n,n,1.d0,eigvectors,size(eigvectors,1), &
+ eigvalues_mat,size(eigvalues_mat,1),0.d0,matrix_tmp1,size(matrix_tmp1,1))
+ ! matrix_tmp2 = W t^-1 sin(t) W^T 
+ call dgemm('N','T',n,n,n,-1.d0,matrix_tmp1,size(matrix_tmp1,1), &
+ eigvectors,size(eigvectors,1),0.d0,matrix_tmp2,size(matrix_tmp2,1))
+ ! exp_X += matrix_tmp2 X 
+ call dgemm('N','N',n,n,n,1.d0,matrix_tmp2,size(matrix_tmp2,1), &
+ X,size(X,1),1.d0,exp_X,size(exp_X,1))
+ 
+end
+
+
+subroutine exp_matrix_taylor(X,n,exp_X,converged)
+ implicit none
+ BEGIN_DOC
+ ! exponential of a general real matrix X using the Taylor expansion of exp(X) 
+ ! 
+ ! returns the logical converged which checks the convergence 
+ END_DOC
+ double precision, intent(in) :: X(n,n)
+ integer, intent(in):: n
+ double precision, intent(out):: exp_X(n,n)
+ logical :: converged
+ double precision :: f
+ integer :: i,iter
+ double precision, allocatable :: Tpotmat(:,:),Tpotmat2(:,:)
+ allocate(Tpotmat(n,n),Tpotmat2(n,n))
+ BEGIN_DOC
+ ! exponential of X using Taylor expansion 
+ END_DOC
+ Tpotmat(:,:)=0.D0
+ exp_X(:,:)  =0.D0
+ do i=1,n
+   Tpotmat(i,i)=1.D0
+   exp_X(i,i)  =1.d0
+ end do
+ iter=0
+ converged=.false.
+ do while (.not.converged)
+   iter+=1
+   f = 1.d0 / dble(iter)
+   Tpotmat2(:,:) = Tpotmat(:,:) * f
+   call dgemm('N','N', n,n,n,1.d0,                   &
+       Tpotmat2, size(Tpotmat2,1),                                  &
+       X, size(X,1), 0.d0,                                    &
+       Tpotmat, size(Tpotmat,1))
+   exp_X(:,:) = exp_X(:,:) + Tpotmat(:,:)
+   
+   converged = ( sum(abs(Tpotmat(:,:))) < 1.d-6).or.(iter>30)
+ end do
+ if(.not.converged)then
+  print*,'Warning !! exp_matrix_taylor did not converge !'
+ endif
+  
+end
+
+subroutine get_A_squared(A,n,A2)
+ implicit none
+ BEGIN_DOC
+! A2 = A A where A is n x n matrix. Use the dgemm routine 
+ END_DOC
+ double precision, intent(in) :: A(n,n)
+ integer, intent(in) :: n
+ double precision, intent(out):: A2(n,n) 
+ call dgemm('N','N',n,n,n,1.d0,A,size(A,1),A,size(A,1),0.d0,A2,size(A2,1))
+end