From fbb946d8f44161bb8de5c752060e82980265717b Mon Sep 17 00:00:00 2001 From: eginer Date: Thu, 15 Feb 2024 16:46:05 +0100 Subject: [PATCH 1/4] removed the systematic save of MOs in casscf --- src/casscf_cipsi/casscf.irp.f | 2 +- src/two_body_rdm/act_2_transition_rdm.irp.f | 4 ++-- 2 files changed, 3 insertions(+), 3 deletions(-) diff --git a/src/casscf_cipsi/casscf.irp.f b/src/casscf_cipsi/casscf.irp.f index addca236..c0cd361d 100644 --- a/src/casscf_cipsi/casscf.irp.f +++ b/src/casscf_cipsi/casscf.irp.f @@ -99,8 +99,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 diff --git a/src/two_body_rdm/act_2_transition_rdm.irp.f b/src/two_body_rdm/act_2_transition_rdm.irp.f index 3d08b084..612213e2 100644 --- a/src/two_body_rdm/act_2_transition_rdm.irp.f +++ b/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'} +! \sum_{\sigma,\sigma'} ! ! WHERE ALL ORBITALS (i,j,k,l) BELONGS TO AN ACTIVE SPACE DEFINED BY "list_act" ! From 22c99a0484eb75ed85c789fa7e39bc965c7fd591 Mon Sep 17 00:00:00 2001 From: eginer Date: Thu, 15 Feb 2024 19:32:15 +0100 Subject: [PATCH 2/4] done some cleaning in the casscf and added a detailed example of Multi state CASSCF --- src/casscf_cipsi/README.rst | 9 ++- src/casscf_cipsi/casscf.irp.f | 16 ++++- src/casscf_cipsi/example_casscf_multistate.sh | 66 +++++++++++++++++++ 3 files changed, 85 insertions(+), 6 deletions(-) create mode 100755 src/casscf_cipsi/example_casscf_multistate.sh diff --git a/src/casscf_cipsi/README.rst b/src/casscf_cipsi/README.rst index f84cde75..75c99de2 100644 --- a/src/casscf_cipsi/README.rst +++ b/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 + diff --git a/src/casscf_cipsi/casscf.irp.f b/src/casscf_cipsi/casscf.irp.f index c0cd361d..d0a26d36 100644 --- a/src/casscf_cipsi/casscf.irp.f +++ b/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 diff --git a/src/casscf_cipsi/example_casscf_multistate.sh b/src/casscf_cipsi/example_casscf_multistate.sh new file mode 100755 index 00000000..368c0440 --- /dev/null +++ b/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 +# 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 +# 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 +# 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} + +# 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 From fa877df399a918750c28a0a262d27823c0cbd3c6 Mon Sep 17 00:00:00 2001 From: eginer Date: Sun, 18 Feb 2024 15:12:39 +0100 Subject: [PATCH 3/4] added exponential of anti-hermitian matrices using the Helgaker's book formulation, and of general matrices using the Taylor expansion. Replaced in casscf_cipsi Umat variable --- src/casscf_cipsi/neworbs.irp.f | 41 +++++----- src/utils/linear_algebra.irp.f | 137 +++++++++++++++++++++++++++++++++ 2 files changed, 158 insertions(+), 20 deletions(-) diff --git a/src/casscf_cipsi/neworbs.irp.f b/src/casscf_cipsi/neworbs.irp.f index a7cebbb2..ca2deebb 100644 --- a/src/casscf_cipsi/neworbs.irp.f +++ b/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 diff --git a/src/utils/linear_algebra.irp.f b/src/utils/linear_algebra.irp.f index 2db47092..175beff3 100644 --- a/src/utils/linear_algebra.irp.f +++ b/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 From ac805f9f016ab5b035557642e19602144d845c6f Mon Sep 17 00:00:00 2001 From: eginer Date: Sun, 18 Feb 2024 15:25:38 +0100 Subject: [PATCH 4/4] added some reference numbers in the example_casscf_multistate.sh --- src/casscf_cipsi/example_casscf_multistate.sh | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/src/casscf_cipsi/example_casscf_multistate.sh b/src/casscf_cipsi/example_casscf_multistate.sh index 368c0440..716c211a 100755 --- a/src/casscf_cipsi/example_casscf_multistate.sh +++ b/src/casscf_cipsi/example_casscf_multistate.sh @@ -10,7 +10,7 @@ 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 +qp run scf # ROHF energy : -149.619992871398 # Freeze the 1s orbitals of the two oxygen qp set_frozen_core @@ -18,7 +18,7 @@ qp set_frozen_core # 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 +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 ! @@ -32,7 +32,7 @@ qp run save_natorb | tee ${EZFIO_FILE}.natorb_3states.out # !!!!! 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 +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 @@ -47,7 +47,7 @@ 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} +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