diff --git a/src/casscf/bielec.irp.f b/src/casscf/bielec.irp.f index 9bb953f8..74351760 100644 --- a/src/casscf/bielec.irp.f +++ b/src/casscf/bielec.irp.f @@ -55,7 +55,6 @@ end do end do - write(6,*) ' provided integrals (PQ|xx) ' END_PROVIDER @@ -116,7 +115,6 @@ BEGIN_PROVIDER [real*8, bielec_PxxQ, (mo_num,n_core_orb+n_act_orb,n_core_orb+n_a end do end do end do - write(6,*) ' provided integrals (Px|xQ) ' END_PROVIDER @@ -146,6 +144,5 @@ BEGIN_PROVIDER [real*8, bielecCI, (n_act_orb,n_act_orb,n_act_orb, mo_num)] end do end do end do - write(6,*) ' provided integrals (tu|xP) ' END_PROVIDER diff --git a/src/casscf/bielec_natorb.irp.f b/src/casscf/bielec_natorb.irp.f index 2f1e43eb..ca1c8e9d 100644 --- a/src/casscf/bielec_natorb.irp.f +++ b/src/casscf/bielec_natorb.irp.f @@ -84,7 +84,6 @@ end do end do end do - write(6,*) ' transformed PQxx' END_PROVIDER @@ -176,7 +175,6 @@ BEGIN_PROVIDER [real*8, bielec_PxxQ_no, (mo_num,n_core_orb+n_act_orb,n_core_orb+ end do end do end do - write(6,*) ' transformed PxxQ ' END_PROVIDER @@ -267,7 +265,6 @@ BEGIN_PROVIDER [real*8, bielecCI_no, (n_act_orb,n_act_orb,n_act_orb, mo_num)] end do end do end do - write(6,*) ' transformed tuvP ' END_PROVIDER diff --git a/src/casscf/casscf.irp.f b/src/casscf/casscf.irp.f index 16c34131..10a3e34a 100644 --- a/src/casscf/casscf.irp.f +++ b/src/casscf/casscf.irp.f @@ -12,20 +12,32 @@ subroutine run implicit none double precision :: energy_old, energy logical :: converged + integer :: iteration converged = .False. energy = 0.d0 -! do while (.not.converged) - N_det = 1 - TOUCH N_det psi_det psi_coef + mo_label = "MCSCF" + iteration = 1 + do while (.not.converged) call run_cipsi - write(6,*) ' total energy = ',eone+etwo+ecore - - call driver_optorb energy_old = energy energy = eone+etwo+ecore - converged = dabs(energy - energy_old) < 1.d-10 -! enddo + + 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_improvement, 'Predicted energy improvement') + + converged = dabs(energy_improvement) < thresh_scf + + mo_coef = NewOrbs + call save_mos + call map_deinit(mo_integrals_map) + N_det = 1 + iteration += 1 + FREE mo_integrals_map mo_two_e_integrals_in_map psi_det psi_coef + SOFT_TOUCH mo_coef N_det + enddo end diff --git a/src/casscf/densities.irp.f b/src/casscf/densities.irp.f index 6e8065e2..8be2db6e 100644 --- a/src/casscf/densities.irp.f +++ b/src/casscf/densities.irp.f @@ -22,7 +22,9 @@ BEGIN_PROVIDER [real*8, D0tu, (n_act_orb,n_act_orb) ] integer :: ierr1,ierr2 real*8 :: cI_mu(N_states) - write(6,*) ' providing density matrices D0 and P0 ' + if (bavard) then + write(6,*) ' providing density matrix D0' + endif D0tu = 0.d0 @@ -90,7 +92,9 @@ BEGIN_PROVIDER [real*8, P0tuvx, (n_act_orb,n_act_orb,n_act_orb,n_act_orb) ] integer(bit_kind), dimension(N_int,2) :: det_mu_ex1, det_mu_ex11, det_mu_ex12 integer(bit_kind), dimension(N_int,2) :: det_mu_ex2, det_mu_ex21, det_mu_ex22 - write(6,*) ' providing density matrices D0 and P0 ' + if (bavard) then + write(6,*) ' providing density matrix P0' + endif P0tuvx = 0.d0 diff --git a/src/casscf/driver_optorb.irp.f b/src/casscf/driver_optorb.irp.f index 591c90c9..2e3e02dc 100644 --- a/src/casscf/driver_optorb.irp.f +++ b/src/casscf/driver_optorb.irp.f @@ -1,32 +1,3 @@ - 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 +subroutine driver_optorb + implicit none +end diff --git a/src/casscf/gradient.irp.f b/src/casscf/gradient.irp.f index 606bf12b..883a4665 100644 --- a/src/casscf/gradient.irp.f +++ b/src/casscf/gradient.irp.f @@ -6,7 +6,6 @@ BEGIN_PROVIDER [ integer, nMonoEx ] 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)] @@ -87,9 +86,11 @@ BEGIN_PROVIDER [real*8, gradvec, (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,*) + if (bavard) then + write(6,*) + write(6,*) ' Norm of the orbital gradient (via <0|EH|0>) : ', norm_grad + write(6,*) + endif END_PROVIDER @@ -118,17 +119,11 @@ subroutine calc_grad_elem(ihole,ipart,res) 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 @@ -176,9 +171,11 @@ BEGIN_PROVIDER [real*8, gradvec2, (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,*) + if (bavard) then + write(6,*) + write(6,*) ' Norm of the orbital gradient (via D, P and integrals): ', norm_grad + write(6,*) + endif END_PROVIDER diff --git a/src/casscf/hessian.irp.f b/src/casscf/hessian.irp.f index 65734a25..e047c5fd 100644 --- a/src/casscf/hessian.irp.f +++ b/src/casscf/hessian.irp.f @@ -14,8 +14,10 @@ BEGIN_PROVIDER [real*8, hessmat, (nMonoEx,nMonoEx)] character*3 :: iexc,jexc real*8 :: res - write(6,*) ' providing Hessian matrix hessmat ' - write(6,*) ' nMonoEx = ',nMonoEx + if (bavard) then + write(6,*) ' providing Hessian matrix hessmat ' + write(6,*) ' nMonoEx = ',nMonoEx + endif do indx=1,nMonoEx do jndx=1,nMonoEx @@ -32,8 +34,6 @@ BEGIN_PROVIDER [real*8, hessmat, (nMonoEx,nMonoEx)] 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 @@ -198,8 +198,10 @@ BEGIN_PROVIDER [real*8, hessmat2, (nMonoEx,nMonoEx)] real*8 :: hessmat_iatb real*8 :: hessmat_taub - write(6,*) ' providing Hessian matrix hessmat2 ' - write(6,*) ' nMonoEx = ',nMonoEx + if (bavard) then + write(6,*) ' providing Hessian matrix hessmat2 ' + write(6,*) ' nMonoEx = ',nMonoEx + endif indx=1 do i=1,n_core_orb @@ -214,7 +216,6 @@ BEGIN_PROVIDER [real*8, hessmat2, (nMonoEx,nMonoEx)] 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 @@ -294,7 +295,6 @@ real*8 function hessmat_itju(i,t,j,u) 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 @@ -340,8 +340,6 @@ real*8 function hessmat_itju(i,t,j,u) end do end do end do - !!! write(6,*) ' direct diff ',i,t,j,u,term,term2 - !!! term=term2 end if else ! it/ju @@ -382,7 +380,6 @@ real*8 function hessmat_itja(i,t,j,a) 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) @@ -416,7 +413,6 @@ real*8 function hessmat_itua(i,t,u,a) 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 @@ -457,7 +453,6 @@ real*8 function hessmat_iajb(i,a,j,b) 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) @@ -495,7 +490,6 @@ real*8 function hessmat_iatb(i,a,t,b) 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) @@ -552,7 +546,6 @@ real*8 function hessmat_taub(t,a,u,b) 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 diff --git a/src/casscf/natorb.irp.f b/src/casscf/natorb.irp.f index 00c9564c..52cd3747 100644 --- a/src/casscf/natorb.irp.f +++ b/src/casscf/natorb.irp.f @@ -14,10 +14,12 @@ occnum(list_act(i))=occ_act(n_act_orb-i+1) end do - write(6,*) ' occupation numbers ' - do i=1,mo_num - write(6,*) i,occnum(i) - end do + if (bavard) then + write(6,*) ' occupation numbers ' + do i=1,mo_num + write(6,*) i,occnum(i) + end do + endif END_PROVIDER @@ -32,14 +34,12 @@ END_PROVIDER call lapack_diag(occ_act,natorbsCI,D0tu,n_act_orb,n_act_orb) - write(6,*) ' found occupation numbers as ' - do i=1,n_act_orb - write(6,*) i,occ_act(i) - end do - if (bavard) then - ! - + write(6,*) ' found occupation numbers as ' + do i=1,n_act_orb + write(6,*) i,occ_act(i) + end do + integer :: nmx real*8 :: xmx do i=1,n_act_orb @@ -152,7 +152,6 @@ BEGIN_PROVIDER [real*8, P0tuvx_no, (n_act_orb,n_act_orb,n_act_orb,n_act_orb)] end do end do end do - write(6,*) ' transformed P0tuvx ' END_PROVIDER @@ -198,7 +197,6 @@ BEGIN_PROVIDER [real*8, one_ints_no, (mo_num,mo_num)] one_ints_no(j,list_act(p))=d(p) end do end do - write(6,*) ' transformed one_ints ' END_PROVIDER @@ -226,148 +224,5 @@ BEGIN_PROVIDER [real*8, NatOrbsFCI, (ao_num,mo_num)] NatOrbsFCI(j,list_act(p))=d(p) end do end do - write(6,*) ' transformed orbitals ' END_PROVIDER - - - - - -subroutine trf_to_natorb() - implicit none - BEGIN_DOC - ! save the diagonal somewhere, in inverse order - ! 4-index-transform the 2-particle density matrix over active orbitals - ! correct the bielectronic integrals - ! correct the monoelectronic integrals - ! put integrals on file, as well orbitals, and the density matrices - ! - END_DOC - integer :: i,j,k,l,t,u,p,q,pp - real*8 :: d(n_act_orb),d1(n_act_orb),d2(n_act_orb) - - ! we recalculate total energies - write(6,*) - write(6,*) ' recalculating energies after the transformation ' - write(6,*) - write(6,*) - real*8 :: e_one_all - real*8 :: e_two_all - integer :: ii - integer :: jj - integer :: t3 - integer :: tt - integer :: u3 - integer :: uu - integer :: v - integer :: v3 - integer :: vv - integer :: x - integer :: x3 - integer :: xx - - e_one_all=0.D0 - e_two_all=0.D0 - do i=1,n_core_orb - ii=list_core(i) - e_one_all+=2.D0*one_ints_no(ii,ii) - do j=1,n_core_orb - jj=list_core(j) - e_two_all+=2.D0*bielec_PQxx_no(ii,ii,j,j)-bielec_PQxx_no(ii,jj,j,i) - end do - do t=1,n_act_orb - tt=list_act(t) - t3=t+n_core_orb - e_two_all += occnum(list_act(t)) * & - (2.d0*bielec_PQxx_no(tt,tt,i,i) - bielec_PQxx_no(tt,ii,i,t3)) - end do - end do - - - - do t=1,n_act_orb - tt=list_act(t) - e_one_all += occnum(list_act(t))*one_ints_no(tt,tt) - do u=1,n_act_orb - uu=list_act(u) - do v=1,n_act_orb - v3=v+n_core_orb - do x=1,n_act_orb - x3=x+n_core_orb - e_two_all +=P0tuvx_no(t,u,v,x)*bielec_PQxx_no(tt,uu,v3,x3) - end do - end do - end do - end do - write(6,*) ' e_one_all = ',e_one_all - write(6,*) ' e_two_all = ',e_two_all - ecore =nuclear_repulsion - ecore_bis=nuclear_repulsion - do i=1,n_core_orb - ii=list_core(i) - ecore +=2.D0*one_ints_no(ii,ii) - ecore_bis+=2.D0*one_ints_no(ii,ii) - do j=1,n_core_orb - jj=list_core(j) - ecore +=2.D0*bielec_PQxx_no(ii,ii,j,j)-bielec_PQxx_no(ii,jj,j,i) - ecore_bis+=2.D0*bielec_PxxQ_no(ii,i,j,jj)-bielec_PxxQ_no(ii,j,j,ii) - end do - end do - eone =0.D0 - eone_bis=0.D0 - etwo =0.D0 - etwo_bis=0.D0 - etwo_ter=0.D0 - do t=1,n_act_orb - tt=list_act(t) - t3=t+n_core_orb - eone += occnum(list_act(t))*one_ints_no(tt,tt) - eone_bis += occnum(list_act(t))*one_ints_no(tt,tt) - do i=1,n_core_orb - ii=list_core(i) - eone += occnum(list_act(t)) * & - (2.D0*bielec_PQxx_no(tt,tt,i,i ) - bielec_PQxx_no(tt,ii,i,t3)) - eone_bis += occnum(list_act(t)) * & - (2.D0*bielec_PxxQ_no(tt,t3,i,ii) - bielec_PxxQ_no(tt,i ,i,tt)) - end do - do u=1,n_act_orb - uu=list_act(u) - u3=u+n_core_orb - 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 - real*8 :: h1,h2,h3 - h1=bielec_PQxx_no(tt,uu,v3,x3) - h2=bielec_PxxQ_no(tt,u3,v3,xx) - h3=bielecCI_no(t,u,v,xx) - etwo +=P0tuvx_no(t,u,v,x)*h1 - etwo_bis+=P0tuvx_no(t,u,v,x)*h2 - etwo_ter+=P0tuvx_no(t,u,v,x)*h3 - if ((abs(h1-h2).gt.1.D-14).or.(abs(h1-h3).gt.1.D-14)) then - write(6,9901) t,u,v,x,h1,h2,h3 - 9901 format('aie: ',4I4,3E20.12) - end if - end do - end do - end do - end do - - write(6,*) ' energy contributions ' - write(6,*) ' core energy = ',ecore,' using PQxx integrals ' - write(6,*) ' core energy (bis) = ',ecore,' using PxxQ integrals ' - write(6,*) ' 1el energy = ',eone ,' using PQxx integrals ' - write(6,*) ' 1el energy (bis) = ',eone ,' using PxxQ integrals ' - write(6,*) ' 2el energy = ',etwo ,' using PQxx integrals ' - write(6,*) ' 2el energy (bis) = ',etwo_bis,' using PxxQ integrals ' - write(6,*) ' 2el energy (ter) = ',etwo_ter,' using tuvP integrals ' - write(6,*) ' ----------------------------------------- ' - write(6,*) ' sum of all = ',eone+etwo+ecore - write(6,*) - SOFT_TOUCH ecore ecore_bis eone eone_bis etwo etwo_bis etwo_ter - -end subroutine trf_to_natorb - diff --git a/src/casscf/neworbs.irp.f b/src/casscf/neworbs.irp.f index fd115880..fd94eb6a 100644 --- a/src/casscf/neworbs.irp.f +++ b/src/casscf/neworbs.irp.f @@ -51,14 +51,16 @@ END_PROVIDER integer :: ierr,matz,i real*8 :: c0 - 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) + if (bavard) then + 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) + endif energy_improvement = SXeigenval(1) integer :: best_vector @@ -73,16 +75,20 @@ END_PROVIDER end if end do - write(6,*) ' SXdiag : eigenvalue for best overlap with ' - write(6,*) ' previous orbitals = ',SXeigenval(best_vector) energy_improvement = SXeigenval(best_vector) + if (bavard) then + write(6,*) ' SXdiag : eigenvalue for best overlap with ' + write(6,*) ' previous orbitals = ',SXeigenval(best_vector) + write(6,*) ' weight of the 1st element ',c0 + endif + 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 diff --git a/src/casscf/tot_en.irp.f b/src/casscf/tot_en.irp.f index 780cd543..ce787232 100644 --- a/src/casscf/tot_en.irp.f +++ b/src/casscf/tot_en.irp.f @@ -42,8 +42,6 @@ end do end do end do - write(6,*) ' e_one_all = ',e_one_all - write(6,*) ' e_two_all = ',e_two_all ecore =nuclear_repulsion ecore_bis=nuclear_repulsion do i=1,n_core_orb @@ -98,24 +96,6 @@ end do end do - write(6,*) ' energy contributions ' - write(6,*) ' core energy = ',ecore,' using PQxx integrals ' - write(6,*) ' core energy (bis) = ',ecore,' using PxxQ integrals ' - write(6,*) ' 1el energy = ',eone ,' using PQxx integrals ' - write(6,*) ' 1el energy (bis) = ',eone ,' using PxxQ integrals ' - write(6,*) ' 2el energy = ',etwo ,' using PQxx integrals ' - write(6,*) ' 2el energy (bis) = ',etwo_bis,' using PxxQ integrals ' - write(6,*) ' 2el energy (ter) = ',etwo_ter,' using tuvP integrals ' - write(6,*) ' ----------------------------------------- ' - write(6,*) ' sum of all = ',eone+etwo+ecore - write(6,*) - write(6,*) ' nuclear (qp) = ',nuclear_repulsion - write(6,*) ' core energy (qp) = ',core_energy - write(6,*) ' 1el energy (qp) = ',psi_energy_h_core(1) - write(6,*) ' 2el energy (qp) = ',psi_energy_two_e(1) - write(6,*) ' nuc + 1 + 2 (qp) = ',nuclear_repulsion+psi_energy_h_core(1)+psi_energy_two_e(1) - write(6,*) ' <0|H|0> (qp) = ',psi_energy_with_nucl_rep(1) - END_PROVIDER