diff --git a/TODO b/TODO index 15bb6d35..abdb618f 100644 --- a/TODO +++ b/TODO @@ -24,6 +24,8 @@ * Example : Simple Hartree-Fock program from scratch * Examples : subroutine example_module +# enleverle psi_det_size for all complicated stuffs with dimension of psi_coef + # Config file for Cray # Documentation de /etc diff --git a/src/casscf/casscf.irp.f b/src/casscf/casscf.irp.f index 6000f734..d83aa271 100644 --- a/src/casscf/casscf.irp.f +++ b/src/casscf/casscf.irp.f @@ -41,7 +41,7 @@ subroutine run mo_occ = occnum call save_mos iteration += 1 - N_det = N_det/2 + N_det = max(N_det/2 ,N_states) psi_det = psi_det_sorted psi_coef = psi_coef_sorted read_wf = .True. diff --git a/src/casscf/get_energy.irp.f b/src/casscf/get_energy.irp.f index 1f572c9e..362da85d 100644 --- a/src/casscf/get_energy.irp.f +++ b/src/casscf/get_energy.irp.f @@ -15,7 +15,13 @@ end subroutine print_grad implicit none - + integer :: i + do i = 1, nMonoEx + if(dabs(gradvec2(i)).gt.1.d-5)then + print*,'' + print*,i,gradvec2(i),excit(:,i) + endif + enddo end subroutine routine_bis diff --git a/src/casscf/grad_old.irp.f b/src/casscf/grad_old.irp.f new file mode 100644 index 00000000..d60a60c8 --- /dev/null +++ b/src/casscf/grad_old.irp.f @@ -0,0 +1,74 @@ + +BEGIN_PROVIDER [real*8, gradvec_old, (nMonoEx)] + BEGIN_DOC + ! calculate the orbital gradient by hand, i.e. for + ! each determinant I we determine the string E_pq |I> (alpha and beta + ! separately) and generate + ! sum_I c_I is then the pq component of the orbital + ! gradient + ! E_pq = a^+_pa_q + a^+_Pa_Q + END_DOC + implicit none + integer :: ii,tt,aa,indx,ihole,ipart,istate + real*8 :: res + + do indx=1,nMonoEx + ihole=excit(1,indx) + ipart=excit(2,indx) + call calc_grad_elem(ihole,ipart,res) + gradvec_old(indx)=res + end do + + real*8 :: norm_grad + norm_grad=0.d0 + do indx=1,nMonoEx + norm_grad+=gradvec_old(indx)*gradvec_old(indx) + end do + norm_grad=sqrt(norm_grad) + if (bavard) then + write(6,*) + write(6,*) ' Norm of the orbital gradient (via <0|EH|0>) : ', norm_grad + write(6,*) + endif + + +END_PROVIDER + +subroutine calc_grad_elem(ihole,ipart,res) + BEGIN_DOC + ! eq 18 of Siegbahn et al, Physica Scripta 1980 + ! we calculate 2 , q=hole, p=particle + END_DOC + implicit none + integer :: ihole,ipart,mu,iii,ispin,ierr,nu,istate + real*8 :: res + integer(bit_kind), allocatable :: det_mu(:,:),det_mu_ex(:,:) + real*8 :: i_H_psi_array(N_states),phase + allocate(det_mu(N_int,2)) + allocate(det_mu_ex(N_int,2)) + + res=0.D0 + + do mu=1,n_det + ! get the string of the determinant + call det_extract(det_mu,mu,N_int) + do ispin=1,2 + ! do the monoexcitation on it + call det_copy(det_mu,det_mu_ex,N_int) + call do_signed_mono_excitation(det_mu,det_mu_ex,nu & + ,ihole,ipart,ispin,phase,ierr) + if (ierr.eq.1) then + call i_H_psi(det_mu_ex,psi_det,psi_coef,N_int & + ,N_det,N_det,N_states,i_H_psi_array) + do istate=1,N_states + res+=i_H_psi_array(istate)*psi_coef(mu,istate)*phase + end do + end if + end do + end do + + ! state-averaged gradient + res*=2.D0/dble(N_states) + +end subroutine calc_grad_elem + diff --git a/src/casscf/gradient.irp.f b/src/casscf/gradient.irp.f index 6bf8b93b..e717e822 100644 --- a/src/casscf/gradient.irp.f +++ b/src/casscf/gradient.irp.f @@ -60,79 +60,6 @@ END_PROVIDER END_PROVIDER -BEGIN_PROVIDER [real*8, gradvec, (nMonoEx)] - BEGIN_DOC - ! calculate the orbital gradient by hand, i.e. for - ! each determinant I we determine the string E_pq |I> (alpha and beta - ! separately) and generate - ! sum_I c_I is then the pq component of the orbital - ! gradient - ! E_pq = a^+_pa_q + a^+_Pa_Q - END_DOC - implicit none - integer :: ii,tt,aa,indx,ihole,ipart,istate - real*8 :: res - - do indx=1,nMonoEx - ihole=excit(1,indx) - ipart=excit(2,indx) - call calc_grad_elem(ihole,ipart,res) - gradvec(indx)=res - end do - - real*8 :: norm_grad - norm_grad=0.d0 - do indx=1,nMonoEx - norm_grad+=gradvec(indx)*gradvec(indx) - end do - norm_grad=sqrt(norm_grad) - if (bavard) then - write(6,*) - write(6,*) ' Norm of the orbital gradient (via <0|EH|0>) : ', norm_grad - write(6,*) - endif - - -END_PROVIDER - -subroutine calc_grad_elem(ihole,ipart,res) - BEGIN_DOC - ! eq 18 of Siegbahn et al, Physica Scripta 1980 - ! we calculate 2 , q=hole, p=particle - END_DOC - implicit none - integer :: ihole,ipart,mu,iii,ispin,ierr,nu,istate - real*8 :: res - integer(bit_kind), allocatable :: det_mu(:,:),det_mu_ex(:,:) - real*8 :: i_H_psi_array(N_states),phase - allocate(det_mu(N_int,2)) - allocate(det_mu_ex(N_int,2)) - - res=0.D0 - - do mu=1,n_det - ! get the string of the determinant - call det_extract(det_mu,mu,N_int) - do ispin=1,2 - ! do the monoexcitation on it - call det_copy(det_mu,det_mu_ex,N_int) - call do_signed_mono_excitation(det_mu,det_mu_ex,nu & - ,ihole,ipart,ispin,phase,ierr) - if (ierr.eq.1) then - call i_H_psi(det_mu_ex,psi_det,psi_coef,N_int & - ,N_det,N_det,N_states,i_H_psi_array) - do istate=1,N_states - res+=i_H_psi_array(istate)*psi_coef(mu,istate)*phase - end do - end if - end do - end do - - ! state-averaged gradient - res*=2.D0/dble(N_states) - -end subroutine calc_grad_elem - BEGIN_PROVIDER [real*8, gradvec2, (nMonoEx)] BEGIN_DOC ! calculate the orbital gradient from density diff --git a/src/casscf/swap_orb.irp.f b/src/casscf/swap_orb.irp.f index 16290782..5d442157 100644 --- a/src/casscf/swap_orb.irp.f +++ b/src/casscf/swap_orb.irp.f @@ -22,13 +22,11 @@ do i = 1, nMonoEx iorder(i) = i vec_tmp(i) = -dabs(SXvector_lowest(i)) - !print*,'vec_tmp(i) = ',i,vec_tmp(i) enddo call dsort(vec_tmp,iorder,nMonoEx) n_max_overlap = 0 do i = 1, nMonoEx if(dabs(vec_tmp(i)).gt.thresh_overlap_switch)then - ! print*,vec_tmp(i),iorder(i) n_max_overlap += 1 max_overlap(n_max_overlap) = iorder(i) endif @@ -107,7 +105,9 @@ endif enddo - print*,'n_orb_swap = ',n_orb_swap + if(n_orb_swap.gt.0)then + print*,'n_orb_swap = ',n_orb_swap + endif do i = 1, n_orb_swap print*,'imono = ',index_orb_swap(i) print*,orb_swap(1,i),'-->',orb_swap(2,i) diff --git a/src/determinants/EZFIO.cfg b/src/determinants/EZFIO.cfg index 1bfe1d10..a8935695 100644 --- a/src/determinants/EZFIO.cfg +++ b/src/determinants/EZFIO.cfg @@ -38,7 +38,7 @@ default: True type: integer doc: Weight used in the calculation of the one-electron density matrix. 0: 1./(c_0^2), 1: 1/N_states, 2: input state-average weight, 3: 1/(Norm_L3(Psi)) interface: ezfio,provider,ocaml -default: 1 +default: 2 [weight_selection] type: integer diff --git a/src/mo_basis/mos.irp.f b/src/mo_basis/mos.irp.f index 610e9a8c..aa04fb01 100644 --- a/src/mo_basis/mos.irp.f +++ b/src/mo_basis/mos.irp.f @@ -91,7 +91,6 @@ BEGIN_PROVIDER [ double precision, mo_coef, (ao_num,mo_num) ] enddo enddo endif - END_PROVIDER BEGIN_PROVIDER [ double precision, mo_coef_in_ao_ortho_basis, (ao_num, mo_num) ] diff --git a/src/two_body_rdm/orb_range_2_rdm.irp.f b/src/two_body_rdm/orb_range_2_rdm.irp.f index b4ec8d87..2bcd04dc 100644 --- a/src/two_body_rdm/orb_range_2_rdm.irp.f +++ b/src/two_body_rdm/orb_range_2_rdm.irp.f @@ -79,7 +79,9 @@ double precision :: wall_0,wall_1 call wall_time(wall_0) print*,'providing the state average TWO-RDM ...' - call orb_range_two_rdm_state_av(state_av_act_two_rdm_spin_trace_mo,n_act_orb,n_act_orb,list_act,list_act_reverse,state_weights,ispin,psi_coef,size(psi_coef,2),size(psi_coef,1)) + print*,'psi_det_size = ',psi_det_size + print*,'N_det = ',N_det + call orb_range_two_rdm_state_av(state_av_act_two_rdm_spin_trace_mo,n_act_orb,n_act_orb,list_act,list_act_reverse,state_weights,ispin,psi_coef,N_states,size(psi_coef,1)) call wall_time(wall_1) print*,'Time to provide the state average TWO-RDM',wall_1 - wall_0 diff --git a/src/two_body_rdm/orb_range_routines.irp.f b/src/two_body_rdm/orb_range_routines.irp.f index a8684185..058ed1c5 100644 --- a/src/two_body_rdm/orb_range_routines.irp.f +++ b/src/two_body_rdm/orb_range_routines.irp.f @@ -147,6 +147,7 @@ subroutine orb_range_two_rdm_state_av_work_$N_int(big_array,dim1,norb,list_orb,l do i=1,maxab idx0(i) = i enddo + ! Prepare the array of all alpha single excitations ! -------------------------------------------------