diff --git a/src/casscf/bielec.irp.f b/src/casscf/bielec.irp.f new file mode 100644 index 00000000..a1ec155d --- /dev/null +++ b/src/casscf/bielec.irp.f @@ -0,0 +1,104 @@ +! -*- F90 -*- + BEGIN_PROVIDER[real*8, bielec_PQxx, (mo_num, mo_num,n_core_orb+n_act_orb,n_core_orb+n_act_orb)] +&BEGIN_PROVIDER[real*8, bielec_PxxQ, (mo_num,n_core_orb+n_act_orb,n_core_orb+n_act_orb, mo_num)] +BEGIN_DOC +! bielec_PQxx : integral (pq|xx) with p,q arbitrary, x core or active +! bielec_PxxQ : integral (px|xq) with p,q arbitrary, x core or active +! indices are unshifted orbital numbers +! all integrals are read from files +END_DOC + implicit none + integer :: i,j,p,q,indx,kk + real*8 :: hhh + logical :: lread + + do i=1,n_core_orb+n_act_orb + do j=1,n_core_orb+n_act_orb + do p=1,mo_num + do q=1,mo_num + bielec_PQxx(p,q,i,j)=0.D0 + bielec_PxxQ(p,i,j,q)=0.D0 + end do + end do + end do + end do + + open(unit=12,form='formatted',status='old',file='bielec_PQxx.tmp') + lread=.true. + indx=0 + do while (lread) + read(12,*,iostat=kk) p,q,i,j,hhh + if (kk.ne.0) then + lread=.false. + else +! stored with p.le.q, i.le.j + bielec_PQxx(p,q,i,j)=hhh + bielec_PQxx(q,p,i,j)=hhh + bielec_PQxx(q,p,j,i)=hhh + bielec_PQxx(p,q,j,i)=hhh + indx+=1 + end if + end do + close(12) + write(6,*) ' read ',indx,' integrals PQxx into core ' + + open(unit=12,form='formatted',status='old',file='bielec_PxxQ.tmp') + lread=.true. + indx=0 + do while (lread) + read(12,*,iostat=kk) p,i,j,q,hhh + if (kk.ne.0) then + lread=.false. + else +! stored with (ip).le.(jq) + bielec_PxxQ(p,i,j,q)=hhh + bielec_PxxQ(q,j,i,p)=hhh + indx+=1 + end if + end do + write(6,*) ' read ',indx,' integrals PxxQ into core ' + close(12) + write(6,*) ' provided integrals (PQ|xx) and (Px|xQ) ' +END_PROVIDER + +BEGIN_PROVIDER[real*8, bielecCI, (n_act_orb,n_act_orb,n_act_orb, mo_num)] +BEGIN_DOC +! bielecCI : integrals (tu|vp) with p arbitrary, tuv active +! index p runs over the whole basis, t,u,v only over the active orbitals +! inegrals read from file +END_DOC + implicit none + integer :: i,j,k,p,t,u,v,kk,indx + real*8 :: hhh + logical :: lread + + write(6,*) ' reading integrals bielecCI ' + + do i=1,n_act_orb + do j=1,n_act_orb + do k=1,n_act_orb + do p=1,mo_num + bielecCI(i,k,j,p)=0.D0 + end do + end do + end do + end do + + open(unit=12,form='formatted',status='old',file='bielecCI.tmp') + lread=.true. + indx=0 + do while (lread) + read(12,*,iostat=kk) i,j,k,p,hhh + if (kk.ne.0) then + lread=.false. + else + bielecCI(i,j,k,p)=hhh + bielecCI(j,i,k,p)=hhh + indx+=1 + end if + end do + write(6,*) ' read ',indx,' integrals (aa|aP) into core ' + close(12) + write(6,*) ' provided integrals (tu|xP) ' +END_PROVIDER + diff --git a/src/casscf/casscf.irp.f b/src/casscf/casscf.irp.f index c08dd032..b55c4c3b 100644 --- a/src/casscf/casscf.irp.f +++ b/src/casscf/casscf.irp.f @@ -10,7 +10,20 @@ end subroutine run implicit none - call run_cipsi - call driver_wdens + double precision :: energy_old, energy + logical :: converged + converged = .False. + + energy = 0.d0 +! do while (.not.converged) + N_det = 1 + TOUCH N_det psi_det psi_coef + call run_cipsi + call driver_wdens + call driver_optorb + energy_old = energy + energy = eone+etwo+ecore + converged = dabs(energy - energy_old) < 1.d-10 +! enddo end diff --git a/src/casscf/driver_optorb.irp.f b/src/casscf/driver_optorb.irp.f new file mode 100644 index 00000000..591c90c9 --- /dev/null +++ b/src/casscf/driver_optorb.irp.f @@ -0,0 +1,32 @@ + subroutine driver_optorb + implicit none + integer :: i,j + + write(6,*) +! write(6,*) ' <0|H|0> (qp) = ',psi_energy_with_nucl_rep(1) + write(6,*) ' energy improvement = ',energy_improvement +! write(6,*) ' new energy = ',psi_energy_with_nucl_rep(1)+energy_improvement + write(6,*) + + write(6,*) + write(6,*) ' creating new orbitals ' + do i=1,mo_num + write(6,*) ' Orbital No ',i + write(6,'(5F14.6)') (NewOrbs(j,i),j=1,mo_num) + write(6,*) + end do + + mo_label = "Natural" + do i=1,mo_num + do j=1,ao_num + mo_coef(j,i)=NewOrbs(j,i) + end do + end do + call save_mos + call map_deinit(mo_integrals_map) + FREE mo_integrals_map mo_coef mo_two_e_integrals_in_map + + write(6,*) + write(6,*) ' ... all done ' + + end diff --git a/src/casscf/gradient.irp.f b/src/casscf/gradient.irp.f new file mode 100644 index 00000000..d35d96ed --- /dev/null +++ b/src/casscf/gradient.irp.f @@ -0,0 +1,251 @@ +! -*- F90 -*- + +use bitmasks ! you need to include the bitmasks_module.f90 features + +BEGIN_PROVIDER [ integer, nMonoEx ] +BEGIN_DOC +! +END_DOC + implicit none + nMonoEx=n_core_orb*n_act_orb+n_core_orb*n_virt_orb+n_act_orb*n_virt_orb + write(6,*) ' nMonoEx = ',nMonoEx +END_PROVIDER + + BEGIN_PROVIDER [integer, excit, (2,nMonoEx)] +&BEGIN_PROVIDER [character*3, excit_class, (nMonoEx)] +BEGIN_DOC +! a list of the orbitals involved in the excitation +END_DOC + + implicit none + integer :: i,t,a,ii,tt,aa,indx + indx=0 + do ii=1,n_core_orb + i=list_core(ii) + do tt=1,n_act_orb + t=list_act(tt) + indx+=1 + excit(1,indx)=i + excit(2,indx)=t + excit_class(indx)='c-a' + end do + end do + + do ii=1,n_core_orb + i=list_core(ii) + do aa=1,n_virt_orb + a=list_virt(aa) + indx+=1 + excit(1,indx)=i + excit(2,indx)=a + excit_class(indx)='c-v' + end do + end do + + do tt=1,n_act_orb + t=list_act(tt) + do aa=1,n_virt_orb + a=list_virt(aa) + indx+=1 + excit(1,indx)=t + excit(2,indx)=a + excit_class(indx)='a-v' + end do + end do + + if (bavard) then + write(6,*) ' Filled the table of the Monoexcitations ' + do indx=1,nMonoEx + write(6,*) ' ex ',indx,' : ',excit(1,indx),' -> ' & + ,excit(2,indx),' ',excit_class(indx) + end do + end if + +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) + write(6,*) + write(6,*) ' Norm of the orbital gradient (via <0|EH|0>) : ', norm_grad + write(6,*) + + +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 +! write(6,*) +! write(6,*) ' mu = ',mu +! call print_det(det_mu,N_int) +! write(6,*) ' generated nu = ',nu,' for excitation ',ihole,' -> ',ipart,' ierr = ',ierr,' phase = ',phase,' ispin = ',ispin +! call print_det(det_mu_ex,N_int) + 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 +! write(6,*) ' contribution = ',i_H_psi_array(1)*psi_coef(mu,1)*phase,res + 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 +! matrices and integrals; Siegbahn et al, Phys Scr 1980 +! eqs 14 a,b,c +END_DOC + implicit none + integer :: i,t,a,indx + real*8 :: gradvec_it,gradvec_ia,gradvec_ta + real*8 :: norm_grad + + indx=0 + do i=1,n_core_orb + do t=1,n_act_orb + indx+=1 + gradvec2(indx)=gradvec_it(i,t) + end do + end do + + do i=1,n_core_orb + do a=1,n_virt_orb + indx+=1 + gradvec2(indx)=gradvec_ia(i,a) + end do + end do + + do t=1,n_act_orb + do a=1,n_virt_orb + indx+=1 + gradvec2(indx)=gradvec_ta(t,a) + end do + end do + + norm_grad=0.d0 + do indx=1,nMonoEx + norm_grad+=gradvec2(indx)*gradvec2(indx) + end do + norm_grad=sqrt(norm_grad) + write(6,*) + write(6,*) ' Norm of the orbital gradient (via D, P and integrals): ', norm_grad + write(6,*) + +END_PROVIDER + + real*8 function gradvec_it(i,t) +BEGIN_DOC +! the orbital gradient core -> active +! we assume natural orbitals +END_DOC + implicit none + integer :: i,t + + integer :: ii,tt,v,vv,x,y + integer :: x3,y3 + + ii=list_core(i) + tt=list_act(t) + gradvec_it=2.D0*(Fipq(tt,ii)+Fapq(tt,ii)) + gradvec_it-=occnum(tt)*Fipq(ii,tt) + do v=1,n_act_orb + vv=list_act(v) + do x=1,n_act_orb + x3=x+n_core_orb + do y=1,n_act_orb + y3=y+n_core_orb + gradvec_it-=2.D0*P0tuvx_no(t,v,x,y)*bielec_PQxx(ii,vv,x3,y3) + end do + end do + end do + gradvec_it*=2.D0 + end function gradvec_it + + real*8 function gradvec_ia(i,a) +BEGIN_DOC +! the orbital gradient core -> virtual +END_DOC + implicit none + integer :: i,a,ii,aa + + ii=list_core(i) + aa=list_virt(a) + gradvec_ia=2.D0*(Fipq(aa,ii)+Fapq(aa,ii)) + gradvec_ia*=2.D0 + + end function gradvec_ia + + real*8 function gradvec_ta(t,a) +BEGIN_DOC +! the orbital gradient active -> virtual +! we assume natural orbitals +END_DOC + implicit none + integer :: t,a,tt,aa,v,vv,x,y + + tt=list_act(t) + aa=list_virt(a) + gradvec_ta=0.D0 + gradvec_ta+=occnum(tt)*Fipq(aa,tt) + do v=1,n_act_orb + do x=1,n_act_orb + do y=1,n_act_orb + gradvec_ta+=2.D0*P0tuvx_no(t,v,x,y)*bielecCI(x,y,v,aa) + end do + end do + end do + gradvec_ta*=2.D0 + + end function gradvec_ta + diff --git a/src/casscf/hessian.irp.f b/src/casscf/hessian.irp.f new file mode 100644 index 00000000..4603d11e --- /dev/null +++ b/src/casscf/hessian.irp.f @@ -0,0 +1,639 @@ +! -*- F90 -*- + +use bitmasks ! you need to include the bitmasks_module.f90 features + +BEGIN_PROVIDER [real*8, hessmat, (nMonoEx,nMonoEx)] +BEGIN_DOC +! calculate the orbital hessian 2 +! + + by hand, +! determinant per determinant, as for the gradient +! +! we assume that we have natural active orbitals +END_DOC + implicit none + integer :: indx,ihole,ipart + integer :: jndx,jhole,jpart + character*3 :: iexc,jexc + real*8 :: res + + write(6,*) ' providing Hessian matrix hessmat ' + write(6,*) ' nMonoEx = ',nMonoEx + + do indx=1,nMonoEx + do jndx=1,nMonoEx + hessmat(indx,jndx)=0.D0 + end do + end do + + do indx=1,nMonoEx + ihole=excit(1,indx) + ipart=excit(2,indx) + iexc=excit_class(indx) + do jndx=indx,nMonoEx + jhole=excit(1,jndx) + jpart=excit(2,jndx) + jexc=excit_class(jndx) + call calc_hess_elem(ihole,ipart,jhole,jpart,res) +! write(6,*) ' Hessian ',ihole,'->',ipart & +! ,' (',iexc,')',jhole,'->',jpart,' (',jexc,')',res + hessmat(indx,jndx)=res + hessmat(jndx,indx)=res + end do + end do + +END_PROVIDER + + subroutine calc_hess_elem(ihole,ipart,jhole,jpart,res) +BEGIN_DOC +! eq 19 of Siegbahn et al, Physica Scripta 1980 +! we calculate 2 +! + + +! average over all states is performed. +! no transition between states. +END_DOC + implicit none + integer :: ihole,ipart,ispin,mu,istate + integer :: jhole,jpart,jspin + integer :: mu_pq, mu_pqrs, mu_rs, mu_rspq, nu_rs,nu + real*8 :: res + integer(bit_kind), allocatable :: det_mu(:,:) + integer(bit_kind), allocatable :: det_nu(:,:) + integer(bit_kind), allocatable :: det_mu_pq(:,:) + integer(bit_kind), allocatable :: det_mu_rs(:,:) + integer(bit_kind), allocatable :: det_nu_rs(:,:) + integer(bit_kind), allocatable :: det_mu_pqrs(:,:) + integer(bit_kind), allocatable :: det_mu_rspq(:,:) + real*8 :: i_H_psi_array(N_states),phase,phase2,phase3 + real*8 :: i_H_j_element + allocate(det_mu(N_int,2)) + allocate(det_nu(N_int,2)) + allocate(det_mu_pq(N_int,2)) + allocate(det_mu_rs(N_int,2)) + allocate(det_nu_rs(N_int,2)) + allocate(det_mu_pqrs(N_int,2)) + allocate(det_mu_rspq(N_int,2)) + integer :: mu_pq_possible + integer :: mu_rs_possible + integer :: nu_rs_possible + integer :: mu_pqrs_possible + integer :: mu_rspq_possible + + res=0.D0 + +! the terms <0|E E H |0> + 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 pq on it + call det_copy(det_mu,det_mu_pq,N_int) + call do_signed_mono_excitation(det_mu,det_mu_pq,mu_pq & + ,ihole,ipart,ispin,phase,mu_pq_possible) + if (mu_pq_possible.eq.1) then +! possible, but not necessarily in the list +! do the second excitation + do jspin=1,2 + call det_copy(det_mu_pq,det_mu_pqrs,N_int) + call do_signed_mono_excitation(det_mu_pq,det_mu_pqrs,mu_pqrs & + ,jhole,jpart,jspin,phase2,mu_pqrs_possible) +! excitation possible + if (mu_pqrs_possible.eq.1) then + call i_H_psi(det_mu_pqrs,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*phase2 + end do + end if +! try the de-excitation with opposite sign + call det_copy(det_mu_pq,det_mu_pqrs,N_int) + call do_signed_mono_excitation(det_mu_pq,det_mu_pqrs,mu_pqrs & + ,jpart,jhole,jspin,phase2,mu_pqrs_possible) + phase2=-phase2 +! excitation possible + if (mu_pqrs_possible.eq.1) then + call i_H_psi(det_mu_pqrs,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*phase2 + end do + end if + end do + end if +! exchange the notion of pq and rs +! do the monoexcitation rs on the initial determinant + call det_copy(det_mu,det_mu_rs,N_int) + call do_signed_mono_excitation(det_mu,det_mu_rs,mu_rs & + ,jhole,jpart,ispin,phase2,mu_rs_possible) + if (mu_rs_possible.eq.1) then +! do the second excitation + do jspin=1,2 + call det_copy(det_mu_rs,det_mu_rspq,N_int) + call do_signed_mono_excitation(det_mu_rs,det_mu_rspq,mu_rspq & + ,ihole,ipart,jspin,phase3,mu_rspq_possible) +! excitation possible (of course, the result is outside the CAS) + if (mu_rspq_possible.eq.1) then + call i_H_psi(det_mu_rspq,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)*phase2*phase3 + end do + end if +! we may try the de-excitation, with opposite sign + call det_copy(det_mu_rs,det_mu_rspq,N_int) + call do_signed_mono_excitation(det_mu_rs,det_mu_rspq,mu_rspq & + ,ipart,ihole,jspin,phase3,mu_rspq_possible) + phase3=-phase3 +! excitation possible (of course, the result is outside the CAS) + if (mu_rspq_possible.eq.1) then + call i_H_psi(det_mu_rspq,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)*phase2*phase3 + end do + end if + end do + end if +! +! the operator E H E, we have to do a double loop over the determinants +! we still have the determinant mu_pq and the phase in memory + if (mu_pq_possible.eq.1) then + do nu=1,N_det + call det_extract(det_nu,nu,N_int) + do jspin=1,2 + call det_copy(det_nu,det_nu_rs,N_int) + call do_signed_mono_excitation(det_nu,det_nu_rs,nu_rs & + ,jhole,jpart,jspin,phase2,nu_rs_possible) +! excitation possible ? + if (nu_rs_possible.eq.1) then + call i_H_j(det_mu_pq,det_nu_rs,N_int,i_H_j_element) + do istate=1,N_states + res+=2.D0*i_H_j_element*psi_coef(mu,istate) & + *psi_coef(nu,istate)*phase*phase2 + end do + end if + end do + end do + end if + end do + end do + +! state-averaged Hessian + res*=1.D0/dble(N_states) + + end subroutine calc_hess_elem + +BEGIN_PROVIDER [real*8, hessmat2, (nMonoEx,nMonoEx)] +BEGIN_DOC +! explicit hessian matrix from density matrices and integrals +! of course, this will be used for a direct Davidson procedure later +! we will not store the matrix in real life +! formulas are broken down as functions for the 6 classes of matrix elements +! +END_DOC + implicit none + integer :: i,j,t,u,a,b,indx,jndx,bstart,ustart + + real*8 :: hessmat_itju + real*8 :: hessmat_itja + real*8 :: hessmat_itua + real*8 :: hessmat_iajb + real*8 :: hessmat_iatb + real*8 :: hessmat_taub + + write(6,*) ' providing Hessian matrix hessmat2 ' + write(6,*) ' nMonoEx = ',nMonoEx + + indx=1 + do i=1,n_core_orb + do t=1,n_act_orb + jndx=indx + do j=i,n_core_orb + if (i.eq.j) then + ustart=t + else + ustart=1 + end if + do u=ustart,n_act_orb + hessmat2(indx,jndx)=hessmat_itju(i,t,j,u) + hessmat2(jndx,indx)=hessmat2(indx,jndx) +! write(6,*) ' result I :',i,t,j,u,indx,jndx,hessmat(indx,jndx),hessmat2(indx,jndx) + jndx+=1 + end do + end do + do j=1,n_core_orb + do a=1,n_virt_orb + hessmat2(indx,jndx)=hessmat_itja(i,t,j,a) + hessmat2(jndx,indx)=hessmat2(indx,jndx) + jndx+=1 + end do + end do + do u=1,n_act_orb + do a=1,n_virt_orb + hessmat2(indx,jndx)=hessmat_itua(i,t,u,a) + hessmat2(jndx,indx)=hessmat2(indx,jndx) + jndx+=1 + end do + end do + indx+=1 + end do + end do + + do i=1,n_core_orb + do a=1,n_virt_orb + jndx=indx + do j=i,n_core_orb + if (i.eq.j) then + bstart=a + else + bstart=1 + end if + do b=bstart,n_virt_orb + hessmat2(indx,jndx)=hessmat_iajb(i,a,j,b) + hessmat2(jndx,indx)=hessmat2(indx,jndx) + jndx+=1 + end do + end do + do t=1,n_act_orb + do b=1,n_virt_orb + hessmat2(indx,jndx)=hessmat_iatb(i,a,t,b) + hessmat2(jndx,indx)=hessmat2(indx,jndx) + jndx+=1 + end do + end do + indx+=1 + end do + end do + + do t=1,n_act_orb + do a=1,n_virt_orb + jndx=indx + do u=t,n_act_orb + if (t.eq.u) then + bstart=a + else + bstart=1 + end if + do b=bstart,n_virt_orb + hessmat2(indx,jndx)=hessmat_taub(t,a,u,b) + hessmat2(jndx,indx)=hessmat2(indx,jndx) + jndx+=1 + end do + end do + indx+=1 + end do + end do + +END_PROVIDER + + real*8 function hessmat_itju(i,t,j,u) +BEGIN_DOC +! the orbital hessian for core->act,core->act +! i, t, j, u are list indices, the corresponding orbitals are ii,tt,jj,uu +! +! we assume natural orbitals +END_DOC + implicit none + integer :: i,t,j,u,ii,tt,uu,v,vv,x,xx,y,jj + real*8 :: term,t2 + +! write(6,*) ' hessmat_itju ',i,t,j,u + ii=list_core(i) + tt=list_act(t) + if (i.eq.j) then + if (t.eq.u) then +! diagonal element + term=occnum(tt)*Fipq(ii,ii)+2.D0*(Fipq(tt,tt)+Fapq(tt,tt)) & + -2.D0*(Fipq(ii,ii)+Fapq(ii,ii)) + term+=2.D0*(3.D0*bielec_pxxq(tt,i,i,tt)-bielec_pqxx(tt,tt,i,i)) + term-=2.D0*occnum(tt)*(3.D0*bielec_pxxq(tt,i,i,tt) & + -bielec_pqxx(tt,tt,i,i)) + term-=occnum(tt)*Fipq(tt,tt) + do v=1,n_act_orb + vv=list_act(v) + do x=1,n_act_orb + xx=list_act(x) + term+=2.D0*(P0tuvx_no(t,t,v,x)*bielec_pqxx(vv,xx,i,i) & + +(P0tuvx_no(t,x,v,t)+P0tuvx_no(t,x,t,v))* & + bielec_pxxq(vv,i,i,xx)) + do y=1,n_act_orb + term-=2.D0*P0tuvx_no(t,v,x,y)*bielecCI(t,v,y,xx) + end do + end do + end do + else +! it/iu, t != u + uu=list_act(u) + term=2.D0*(Fipq(tt,uu)+Fapq(tt,uu)) + term+=2.D0*(4.D0*bielec_PxxQ(tt,i,j,uu)-bielec_PxxQ(uu,i,j,tt) & + -bielec_PQxx(tt,uu,i,j)) + term-=occnum(tt)*Fipq(uu,tt) + term-=(occnum(tt)+occnum(uu)) & + *(3.D0*bielec_PxxQ(tt,i,i,uu)-bielec_PQxx(uu,tt,i,i)) + do v=1,n_act_orb + vv=list_act(v) +! term-=D0tu(u,v)*Fipq(tt,vv) ! published, but inverting t and u seems more correct + do x=1,n_act_orb + xx=list_act(x) + term+=2.D0*(P0tuvx_no(u,t,v,x)*bielec_pqxx(vv,xx,i,i) & + +(P0tuvx_no(u,x,v,t)+P0tuvx_no(u,x,t,v)) & + *bielec_pxxq(vv,i,i,xx)) + do y=1,n_act_orb + term-=2.D0*P0tuvx_no(t,v,x,y)*bielecCI(u,v,y,xx) + end do + end do + end do +!!! write(6,*) ' direct diff ',i,t,j,u,term,term2 +!!! term=term2 + end if + else +! it/ju + jj=list_core(j) + uu=list_act(u) + if (t.eq.u) then + term=occnum(tt)*Fipq(ii,jj) + term-=2.D0*(Fipq(ii,jj)+Fapq(ii,jj)) + else + term=0.D0 + end if + term+=2.D0*(4.D0*bielec_PxxQ(tt,i,j,uu)-bielec_PxxQ(uu,i,j,tt) & + -bielec_PQxx(tt,uu,i,j)) + term-=(occnum(tt)+occnum(uu))* & + (4.D0*bielec_PxxQ(tt,i,j,uu)-bielec_PxxQ(uu,i,j,tt) & + -bielec_PQxx(uu,tt,i,j)) + do v=1,n_act_orb + vv=list_act(v) + do x=1,n_act_orb + xx=list_act(x) + term+=2.D0*(P0tuvx_no(u,t,v,x)*bielec_pqxx(vv,xx,i,j) & + +(P0tuvx_no(u,x,v,t)+P0tuvx_no(u,x,t,v)) & + *bielec_pxxq(vv,i,j,xx)) + end do + end do + end if + + term*=2.D0 + hessmat_itju=term + + end function hessmat_itju + + real*8 function hessmat_itja(i,t,j,a) +BEGIN_DOC +! the orbital hessian for core->act,core->virt +END_DOC + implicit none + integer :: i,t,j,a,ii,tt,jj,aa,v,vv,x,y + real*8 :: term + +! write(6,*) ' hessmat_itja ',i,t,j,a +! it/ja + ii=list_core(i) + tt=list_act(t) + jj=list_core(j) + aa=list_virt(a) + term=2.D0*(4.D0*bielec_pxxq(aa,j,i,tt) & + -bielec_pqxx(aa,tt,i,j) -bielec_pxxq(aa,i,j,tt)) + term-=occnum(tt)*(4.D0*bielec_pxxq(aa,j,i,tt) & + -bielec_pqxx(aa,tt,i,j) -bielec_pxxq(aa,i,j,tt)) + if (i.eq.j) then + term+=2.D0*(Fipq(aa,tt)+Fapq(aa,tt)) + term-=0.5D0*occnum(tt)*Fipq(aa,tt) + do v=1,n_act_orb + do x=1,n_act_orb + do y=1,n_act_orb + term-=P0tuvx_no(t,v,x,y)*bielecCI(x,y,v,aa) + end do + end do + end do + end if + term*=2.D0 + hessmat_itja=term + + end function hessmat_itja + + real*8 function hessmat_itua(i,t,u,a) +BEGIN_DOC +! the orbital hessian for core->act,act->virt +END_DOC + implicit none + integer :: i,t,u,a,ii,tt,uu,aa,v,vv,x,xx,u3,t3,v3 + real*8 :: term + +! write(6,*) ' hessmat_itua ',i,t,u,a + ii=list_core(i) + tt=list_act(t) + t3=t+n_core_orb + uu=list_act(u) + u3=u+n_core_orb + aa=list_virt(a) + if (t.eq.u) then + term=-occnum(tt)*Fipq(aa,ii) + else + term=0.D0 + end if + term-=occnum(uu)*(bielec_pqxx(aa,ii,t3,u3)-4.D0*bielec_pqxx(aa,uu,t3,i) & + +bielec_pxxq(aa,t3,u3,ii)) + do v=1,n_act_orb + vv=list_act(v) + v3=v+n_core_orb + do x=1,n_act_orb +integer :: x3 + xx=list_act(x) + x3=x+n_core_orb + term-=2.D0*(P0tuvx_no(t,u,v,x)*bielec_pqxx(aa,ii,v3,x3) & + +(P0tuvx_no(t,v,u,x)+P0tuvx_no(t,v,x,u)) & + *bielec_pqxx(aa,xx,v3,i)) + end do + end do + if (t.eq.u) then + term+=Fipq(aa,ii)+Fapq(aa,ii) + end if + term*=2.D0 + hessmat_itua=term + + end function hessmat_itua + + real*8 function hessmat_iajb(i,a,j,b) +BEGIN_DOC +! the orbital hessian for core->virt,core->virt +END_DOC + implicit none + integer :: i,a,j,b,ii,aa,jj,bb + real*8 :: term +! write(6,*) ' hessmat_iajb ',i,a,j,b + + ii=list_core(i) + aa=list_virt(a) + if (i.eq.j) then + if (a.eq.b) then +! ia/ia + term=2.D0*(Fipq(aa,aa)+Fapq(aa,aa)-Fipq(ii,ii)-Fapq(ii,ii)) + term+=2.D0*(3.D0*bielec_pxxq(aa,i,i,aa)-bielec_pqxx(aa,aa,i,i)) + else + bb=list_virt(b) +! ia/ib + term=2.D0*(Fipq(aa,bb)+Fapq(aa,bb)) + term+=2.D0*(3.D0*bielec_pxxq(aa,i,i,bb)-bielec_pqxx(aa,bb,i,i)) + end if + else +! ia/jb + jj=list_core(j) + bb=list_virt(b) + term=2.D0*(4.D0*bielec_pxxq(aa,i,j,bb)-bielec_pqxx(aa,bb,i,j) & + -bielec_pxxq(aa,j,i,bb)) + if (a.eq.b) then + term-=2.D0*(Fipq(ii,jj)+Fapq(ii,jj)) + end if + end if + term*=2.D0 + hessmat_iajb=term + + end function hessmat_iajb + + real*8 function hessmat_iatb(i,a,t,b) +BEGIN_DOC +! the orbital hessian for core->virt,act->virt +END_DOC + implicit none + integer :: i,a,t,b,ii,aa,tt,bb,v,vv,x,y,v3,t3 + real*8 :: term + +! write(6,*) ' hessmat_iatb ',i,a,t,b + ii=list_core(i) + aa=list_virt(a) + tt=list_act(t) + bb=list_virt(b) + t3=t+n_core_orb + term=occnum(tt)*(4.D0*bielec_pxxq(aa,i,t3,bb)-bielec_pxxq(aa,t3,i,bb) & + -bielec_pqxx(aa,bb,i,t3)) + if (a.eq.b) then + term-=Fipq(tt,ii)+Fapq(tt,ii) + term-=0.5D0*occnum(tt)*Fipq(tt,ii) + do v=1,n_act_orb + do x=1,n_act_orb + do y=1,n_act_orb + term-=P0tuvx_no(t,v,x,y)*bielecCI(x,y,v,ii) + end do + end do + end do + end if + term*=2.D0 + hessmat_iatb=term + + end function hessmat_iatb + + real*8 function hessmat_taub(t,a,u,b) +BEGIN_DOC +! the orbital hessian for act->virt,act->virt +END_DOC + implicit none + integer :: t,a,u,b,tt,aa,uu,bb,v,vv,x,xx,y + integer :: v3,x3 + real*8 :: term,t1,t2,t3 + + tt=list_act(t) + aa=list_virt(a) + if (t.eq.u) then + if (a.eq.b) then +! ta/ta + t1=occnum(tt)*Fipq(aa,aa) + t2=0.D0 + t3=0.D0 + t1-=occnum(tt)*Fipq(tt,tt) + do v=1,n_act_orb + vv=list_act(v) + v3=v+n_core_orb + do x=1,n_act_orb + xx=list_act(x) + x3=x+n_core_orb + t2+=2.D0*(P0tuvx_no(t,t,v,x)*bielec_pqxx(aa,aa,v3,x3) & + +(P0tuvx_no(t,x,v,t)+P0tuvx_no(t,x,t,v))* & + bielec_pxxq(aa,x3,v3,aa)) + do y=1,n_act_orb + t3-=2.D0*P0tuvx_no(t,v,x,y)*bielecCI(t,v,y,xx) + end do + end do + end do + term=t1+t2+t3 +! write(6,*) ' Hess taub ',t,a,t1,t2,t3 + else + bb=list_virt(b) +! ta/tb b/=a + term=occnum(tt)*Fipq(aa,bb) + do v=1,n_act_orb + vv=list_act(v) + v3=v+n_core_orb + do x=1,n_act_orb + xx=list_act(x) + x3=x+n_core_orb + term+=2.D0*(P0tuvx_no(t,t,v,x)*bielec_pqxx(aa,bb,v3,x3) & + +(P0tuvx_no(t,x,v,t)+P0tuvx_no(t,x,t,v)) & + *bielec_pxxq(aa,x3,v3,bb)) + end do + end do + end if + else +! ta/ub t/=u + uu=list_act(u) + bb=list_virt(b) + term=0.D0 + do v=1,n_act_orb + vv=list_act(v) + v3=v+n_core_orb + do x=1,n_act_orb + xx=list_act(x) + x3=x+n_core_orb + term+=2.D0*(P0tuvx_no(t,u,v,x)*bielec_pqxx(aa,bb,v3,x3) & + +(P0tuvx_no(t,x,v,u)+P0tuvx_no(t,x,u,v)) & + *bielec_pxxq(aa,x3,v3,bb)) + end do + end do + if (a.eq.b) then + term-=0.5D0*(occnum(tt)*Fipq(uu,tt)+occnum(uu)*Fipq(tt,uu)) + do v=1,n_act_orb + do x=1,n_act_orb + do y=1,n_act_orb + term-=P0tuvx_no(t,v,x,y)*bielecCI(x,y,v,uu) + term-=P0tuvx_no(u,v,x,y)*bielecCI(x,y,v,tt) + end do + end do + end do + end if + + end if + + term*=2.D0 + hessmat_taub=term + + end function hessmat_taub + +BEGIN_PROVIDER [real*8, hessdiag, (nMonoEx)] +BEGIN_DOC +! the diagonal of the Hessian, needed for the Davidson procedure +END_DOC + implicit none + integer :: i,t,a,indx + real*8 :: hessmat_itju,hessmat_iajb,hessmat_taub + + indx=0 + do i=1,n_core_orb + do t=1,n_act_orb + indx+=1 + hessdiag(indx)=hessmat_itju(i,t,i,t) + end do + end do + + do i=1,n_core_orb + do a=1,n_virt_orb + indx+=1 + hessdiag(indx)=hessmat_iajb(i,a,i,a) + end do + end do + + do t=1,n_act_orb + do a=1,n_virt_orb + indx+=1 + hessdiag(indx)=hessmat_taub(t,a,t,a) + end do + end do + +END_PROVIDER diff --git a/src/casscf/mcscf_fock.irp.f b/src/casscf/mcscf_fock.irp.f new file mode 100644 index 00000000..301b1418 --- /dev/null +++ b/src/casscf/mcscf_fock.irp.f @@ -0,0 +1,67 @@ +! -*- F90 -*- + BEGIN_PROVIDER [real*8, Fipq, (mo_num,mo_num) ] +&BEGIN_PROVIDER [real*8, Fapq, (mo_num,mo_num) ] +BEGIN_DOC +! the inactive and the active Fock matrices, in molecular +! orbitals +! we create them in MOs, quite expensive +! +! for an implementation in AOs we need first the natural orbitals +! for forming an active density matrix in AOs +! +END_DOC + implicit none + double precision, allocatable :: integrals_array1(:,:) + double precision, allocatable :: integrals_array2(:,:) + integer :: p,q,k,kk,t,tt,u,uu + allocate(integrals_array1(mo_num,mo_num)) + allocate(integrals_array2(mo_num,mo_num)) + + do p=1,mo_num + do q=1,mo_num + Fipq(p,q)=one_ints(p,q) + Fapq(p,q)=0.D0 + end do + end do + +! the inactive Fock matrix + do k=1,n_core_orb + kk=list_core(k) + do p=1,mo_num + do q=1,mo_num + Fipq(p,q)+=2.D0*bielec_pqxx(p,q,k,k) -bielec_pxxq(p,k,k,q) + end do + end do + end do + +! the active Fock matrix, D0tu is diagonal + do t=1,n_act_orb + tt=list_act(t) + do p=1,mo_num + do q=1,mo_num + Fapq(p,q)+=occnum(tt) & + *(bielec_pqxx(p,q,tt,tt)-0.5D0*bielec_pxxq(p,tt,tt,q)) + end do + end do + end do + +if (bavard) then +integer :: i + write(6,*) + write(6,*) ' the effective Fock matrix over MOs' + write(6,*) + + write(6,*) + write(6,*) ' the diagonal of the inactive effective Fock matrix ' + write(6,'(5(i3,F12.5))') (i,Fipq(i,i),i=1,mo_num) + write(6,*) + write(6,*) + write(6,*) ' the diagonal of the active Fock matrix ' + write(6,'(5(i3,F12.5))') (i,Fapq(i,i),i=1,mo_num) + write(6,*) +end if + + +END_PROVIDER + + diff --git a/src/casscf/natorb_casscf.irp.f b/src/casscf/natorb_casscf.irp.f new file mode 100644 index 00000000..0a818a34 --- /dev/null +++ b/src/casscf/natorb_casscf.irp.f @@ -0,0 +1,65 @@ +! -*- F90 -*- +BEGIN_PROVIDER [real*8, occnum, (mo_num)] + implicit none + integer :: i,kk,j + logical :: lread + real*8 :: rdum + do i=1,mo_num + occnum(i)=0.D0 + end do + do i=1,n_core_orb + occnum(list_core(i))=2.D0 + end do + + open(unit=12,file='D0tu.dat',form='formatted',status='old') + lread=.true. + do while (lread) + read(12,*,iostat=kk) i,j,rdum + if (kk.ne.0) then + lread=.false. + else + if (i.eq.j) then + occnum(list_act(i))=rdum + else + write(6,*) ' WARNING - no natural orbitals !' + write(6,*) i,j,rdum + end if + end if + end do + close(12) + write(6,*) ' read occupation numbers ' + do i=1,mo_num + write(6,*) i,occnum(i) + end do + +END_PROVIDER + +BEGIN_PROVIDER [real*8, P0tuvx_no, (n_act_orb,n_act_orb,n_act_orb,n_act_orb)] + implicit none + integer :: i,j,k,l,kk + real*8 :: rdum + logical :: lread + + do i=1,n_act_orb + do j=1,n_act_orb + do k=1,n_act_orb + do l=1,n_act_orb + P0tuvx_no(l,k,j,i)=0.D0 + end do + end do + end do + end do + + open(unit=12,file='P0tuvx.dat',form='formatted',status='old') + lread=.true. + do while (lread) + read(12,*,iostat=kk) i,j,k,l,rdum + if (kk.ne.0) then + lread=.false. + else + P0tuvx_no(i,j,k,l)=rdum + end if + end do + close(12) + write(6,*) ' read the 2-particle density matrix ' +END_PROVIDER diff --git a/src/casscf/neworbs.irp.f b/src/casscf/neworbs.irp.f new file mode 100644 index 00000000..6d63a86e --- /dev/null +++ b/src/casscf/neworbs.irp.f @@ -0,0 +1,222 @@ +! -*- F90 -*- +BEGIN_PROVIDER [real*8, SXmatrix, (nMonoEx+1,nMonoEx+1)] + implicit none + integer :: i,j + do i=1,nMonoEx+1 + do j=1,nMonoEx+1 + SXmatrix(i,j)=0.D0 + end do + end do + + do i=1,nMonoEx + SXmatrix(1,i+1)=gradvec2(i) + SXmatrix(1+i,1)=gradvec2(i) + end do + + do i=1,nMonoEx + do j=1,nMonoEx + SXmatrix(i+1,j+1)=hessmat2(i,j) + SXmatrix(j+1,i+1)=hessmat2(i,j) + end do + end do + + if (bavard) then + do i=2,nMonoEx+1 + write(6,*) ' diagonal of the Hessian : ',i,hessmat2(i,i) + end do + end if + + +END_PROVIDER + + BEGIN_PROVIDER [real*8, SXeigenvec, (nMonoEx+1,nMonoEx+1)] +&BEGIN_PROVIDER [real*8, SXeigenval, (nMonoEx+1)] + END_PROVIDER + + BEGIN_PROVIDER [real*8, SXvector, (nMonoEx+1)] +&BEGIN_PROVIDER [real*8, energy_improvement] + implicit none + integer :: ierr,matz,i + real*8 :: c0 + + call lapack_diag(SXeigenval,SXeigenvec,SXmatrix,nMonoEx+1,nMonoEx+1) + write(6,*) ' SXdiag : lowest 5 eigenvalues ' + write(6,*) ' 1 - ',SXeigenval(1),SXeigenvec(1,1) + write(6,*) ' 2 - ',SXeigenval(2),SXeigenvec(1,2) + write(6,*) ' 3 - ',SXeigenval(3),SXeigenvec(1,3) + write(6,*) ' 4 - ',SXeigenval(4),SXeigenvec(1,4) + write(6,*) ' 5 - ',SXeigenval(5),SXeigenvec(1,5) + write(6,*) + write(6,*) ' SXdiag : lowest eigenvalue = ',SXeigenval(1) + energy_improvement = SXeigenval(1) + +integer :: best_vector +real*8 :: best_overlap + best_overlap=0.D0 + do i=1,nMonoEx+1 + if (SXeigenval(i).lt.0.D0) then + if (abs(SXeigenvec(1,i)).gt.best_overlap) then + best_overlap=abs(SXeigenvec(1,i)) + best_vector=i + end if + end if + end do + + write(6,*) ' SXdiag : eigenvalue for best overlap with ' + write(6,*) ' previous orbitals = ',SXeigenval(best_vector) + energy_improvement = SXeigenval(best_vector) + + c0=SXeigenvec(1,best_vector) + write(6,*) ' weight of the 1st element ',c0 + do i=1,nMonoEx+1 + SXvector(i)=SXeigenvec(i,best_vector)/c0 +! write(6,*) ' component No ',i,' : ',SXvector(i) + end do + +END_PROVIDER + + +BEGIN_PROVIDER [real*8, NewOrbs, (ao_num,mo_num) ] + implicit none + integer :: i,j,ialph + +! form the exponential of the Orbital rotations + call get_orbrotmat +! form the new orbitals + do i=1,ao_num + do j=1,mo_num + NewOrbs(i,j)=0.D0 + end do + end do + + do ialph=1,ao_num + do i=1,mo_num + wrkline(i)=mo_coef(ialph,i) + end do + do i=1,mo_num + do j=1,mo_num + NewOrbs(ialph,i)+=Umat(i,j)*wrkline(j) + end do + end do + end do + +END_PROVIDER + + BEGIN_PROVIDER [real*8, Tpotmat, (mo_num,mo_num) ] +&BEGIN_PROVIDER [real*8, Umat, (mo_num,mo_num) ] +&BEGIN_PROVIDER [real*8, wrkline, (mo_num) ] +&BEGIN_PROVIDER [real*8, Tmat, (mo_num,mo_num) ] +END_PROVIDER + + subroutine get_orbrotmat + implicit none + integer :: i,j,indx,k,iter,t,a,ii,tt,aa + real*8 :: sum + logical :: converged + + +! the orbital rotation matrix T + do i=1,mo_num + do j=1,mo_num + Tmat(i,j)=0.D0 + Umat(i,j)=0.D0 + Tpotmat(i,j)=0.D0 + end do + Tpotmat(i,i)=1.D0 + end do + + indx=1 + do i=1,n_core_orb + ii=list_core(i) + do t=1,n_act_orb + tt=list_act(t) + indx+=1 + Tmat(ii,tt)= SXvector(indx) + Tmat(tt,ii)=-SXvector(indx) + end do + end do + do i=1,n_core_orb + ii=list_core(i) + do a=1,n_virt_orb + aa=list_virt(a) + indx+=1 + Tmat(ii,aa)= SXvector(indx) + Tmat(aa,ii)=-SXvector(indx) + end do + end do + do t=1,n_act_orb + tt=list_act(t) + do a=1,n_virt_orb + aa=list_virt(a) + indx+=1 + Tmat(tt,aa)= SXvector(indx) + Tmat(aa,tt)=-SXvector(indx) + end do + end do + + write(6,*) ' the T matrix ' + do indx=1,nMonoEx + i=excit(1,indx) + j=excit(2,indx) +! if (abs(Tmat(i,j)).gt.1.D0) then +! write(6,*) ' setting matrix element ',i,j,' of ',Tmat(i,j),' to ' & +! , sign(1.D0,Tmat(i,j)) +! Tmat(i,j)=sign(1.D0,Tmat(i,j)) +! Tmat(j,i)=-Tmat(i,j) +! end if + if (abs(Tmat(i,j)).gt.1.D-9) write(6,9901) i,j,excit_class(indx),Tmat(i,j) + 9901 format(' ',i4,' -> ',i4,' (',A3,') : ',E14.6) + end do + + write(6,*) + write(6,*) ' forming the matrix exponential ' + write(6,*) +! form the exponential + iter=0 + converged=.false. + do while (.not.converged) + iter+=1 +! add the next term + do i=1,mo_num + do j=1,mo_num + Umat(i,j)+=Tpotmat(i,j) + end do + end do +! next power of T, we multiply Tpotmat with Tmat/iter + do i=1,mo_num + do j=1,mo_num + wrkline(j)=Tpotmat(i,j)/dble(iter) + Tpotmat(i,j)=0.D0 + end do + do j=1,mo_num + do k=1,mo_num + Tpotmat(i,j)+=wrkline(k)*Tmat(k,j) + end do + end do + end do +! Convergence test + sum=0.D0 + do i=1,mo_num + do j=1,mo_num + sum+=abs(Tpotmat(i,j)) + end do + end do + write(6,*) ' Iteration No ',iter,' Sum = ',sum + if (sum.lt.1.D-6) then + converged=.true. + end if + if (iter.ge.NItExpMax) then + stop ' no convergence ' + end if + end do + write(6,*) + write(6,*) ' Converged ! ' + write(6,*) + + end subroutine get_orbrotmat + +BEGIN_PROVIDER [integer, NItExpMax] + NItExpMax=100 +END_PROVIDER + + diff --git a/src/casscf/one_ints.irp.f b/src/casscf/one_ints.irp.f new file mode 100644 index 00000000..a802f644 --- /dev/null +++ b/src/casscf/one_ints.irp.f @@ -0,0 +1,26 @@ +! -*- F90 -*- +BEGIN_PROVIDER [real*8, one_ints, (mo_num,mo_num)] + implicit none + integer :: i,j,kk + logical :: lread + real*8 :: rdum + do i=1,mo_num + do j=1,mo_num + one_ints(i,j)=0.D0 + end do + end do + open(unit=12,file='onetrf.tmp',status='old',form='formatted') + lread=.true. + do while (lread) + read(12,*,iostat=kk) i,j,rdum + if (kk.ne.0) then + lread=.false. + else + one_ints(i,j)=rdum + one_ints(j,i)=rdum + end if + end do + close(12) + write(6,*) ' read MCSCF natural one-electron integrals ' +END_PROVIDER +