From 92e44f53bae20a9995f968c16017ca402c7c5962 Mon Sep 17 00:00:00 2001 From: Anthony Scemama Date: Thu, 27 Jun 2019 23:06:35 +0200 Subject: [PATCH 1/5] Fixed small bugs --- src/casscf/casscf.irp.f | 5 +++-- src/casscf/neworbs.irp.f | 6 +++--- 2 files changed, 6 insertions(+), 5 deletions(-) diff --git a/src/casscf/casscf.irp.f b/src/casscf/casscf.irp.f index 10a3e34a..8aaa1925 100644 --- a/src/casscf/casscf.irp.f +++ b/src/casscf/casscf.irp.f @@ -3,8 +3,8 @@ program casscf BEGIN_DOC ! TODO : Put the documentation of the program here END_DOC - no_vvvv_integrals = .True. - SOFT_TOUCH no_vvvv_integrals +! no_vvvv_integrals = .True. +! SOFT_TOUCH no_vvvv_integrals call run end @@ -13,6 +13,7 @@ subroutine run double precision :: energy_old, energy logical :: converged integer :: iteration + PROVIDE mo_two_e_integrals_in_map converged = .False. energy = 0.d0 diff --git a/src/casscf/neworbs.irp.f b/src/casscf/neworbs.irp.f index fd94eb6a..f4319485 100644 --- a/src/casscf/neworbs.irp.f +++ b/src/casscf/neworbs.irp.f @@ -25,7 +25,7 @@ BEGIN_PROVIDER [real*8, SXmatrix, (nMonoEx+1,nMonoEx+1)] end do if (bavard) then - do i=2,nMonoEx+1 + do i=2,nMonoEx write(6,*) ' diagonal of the Hessian : ',i,hessmat2(i,i) end do end if @@ -77,14 +77,14 @@ END_PROVIDER energy_improvement = SXeigenval(best_vector) + c0=SXeigenvec(1,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) - do i=1,nMonoEx+1 SXvector(i)=SXeigenvec(i,best_vector)/c0 end do From 82bbf95fead74f927ba1c15779b958623dcfc3ef Mon Sep 17 00:00:00 2001 From: Anthony Scemama Date: Thu, 27 Jun 2019 23:46:30 +0200 Subject: [PATCH 2/5] Fixed small bugs --- src/bitmask/core_inact_act_virt.irp.f | 4 ++++ src/casscf/casscf.irp.f | 9 ++++----- src/cipsi/pt2_stoch_routines.irp.f | 13 +++++++++++-- src/generators_cas/generators.irp.f | 1 + 4 files changed, 20 insertions(+), 7 deletions(-) diff --git a/src/bitmask/core_inact_act_virt.irp.f b/src/bitmask/core_inact_act_virt.irp.f index f830da4e..177c3df5 100644 --- a/src/bitmask/core_inact_act_virt.irp.f +++ b/src/bitmask/core_inact_act_virt.irp.f @@ -141,6 +141,10 @@ END_PROVIDER n_act_orb_tmp = 0 n_virt_orb_tmp = 0 n_del_orb_tmp = 0 + core_bitmask = 0_bit_kind + inact_bitmask = 0_bit_kind + act_bitmask = 0_bit_kind + virt_bitmask = 0_bit_kind do i = 1, mo_num if(mo_class(i) == 'Core')then n_core_orb_tmp += 1 diff --git a/src/casscf/casscf.irp.f b/src/casscf/casscf.irp.f index 8aaa1925..54bf35ee 100644 --- a/src/casscf/casscf.irp.f +++ b/src/casscf/casscf.irp.f @@ -3,8 +3,8 @@ program casscf BEGIN_DOC ! TODO : Put the documentation of the program here END_DOC -! no_vvvv_integrals = .True. -! SOFT_TOUCH no_vvvv_integrals + no_vvvv_integrals = .True. + SOFT_TOUCH no_vvvv_integrals call run end @@ -13,15 +13,13 @@ subroutine run double precision :: energy_old, energy logical :: converged integer :: iteration - PROVIDE mo_two_e_integrals_in_map converged = .False. energy = 0.d0 mo_label = "MCSCF" iteration = 1 do while (.not.converged) - call run_cipsi - + call run_stochastic_cipsi energy_old = energy energy = eone+etwo+ecore @@ -39,6 +37,7 @@ subroutine run 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/cipsi/pt2_stoch_routines.irp.f b/src/cipsi/pt2_stoch_routines.irp.f index 9f891320..7825d24c 100644 --- a/src/cipsi/pt2_stoch_routines.irp.f +++ b/src/cipsi/pt2_stoch_routines.irp.f @@ -135,7 +135,7 @@ subroutine ZMQ_pt2(E, pt2,relative_error, error, variance, norm, N_in) PROVIDE psi_occ_pattern_hii det_to_occ_pattern endif - if (N_det < max(4,N_states)) then + if (N_det <= max(4,N_states)) then pt2=0.d0 variance=0.d0 norm=0.d0 @@ -719,6 +719,15 @@ END_PROVIDER double precision :: rss double precision, external :: memory_of_double, memory_of_int + if (N_det_generators == 1) then + pt2_w = 1.d0 + pt2_cw = 1.d0 + pt2_W_T = 1.d0 + pt2_u_0 = 1.d0 + pt2_n_0 = 1 + return + endif + rss = memory_of_double(2*N_det_generators+1) call check_mem(rss,irp_here) @@ -754,7 +763,7 @@ END_PROVIDER end if pt2_n_0(1) += 1 if(N_det_generators - pt2_n_0(1) < pt2_minDetInFirstTeeth * pt2_N_teeth) then - stop "teeth building failed" + print *, "teeth building failed" end if end do !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! diff --git a/src/generators_cas/generators.irp.f b/src/generators_cas/generators.irp.f index c22eab51..b2f58202 100644 --- a/src/generators_cas/generators.irp.f +++ b/src/generators_cas/generators.irp.f @@ -55,6 +55,7 @@ END_PROVIDER nongen(inongen) = i endif enddo + ASSERT (m == N_det_generators) psi_det_sorted_gen(:,:,:N_det_generators) = psi_det_generators(:,:,:N_det_generators) psi_coef_sorted_gen(:N_det_generators, :) = psi_coef_generators(:N_det_generators, :) From ae3a4929b6443f13512e2c1b28642812914bc84a Mon Sep 17 00:00:00 2001 From: Anthony Scemama Date: Thu, 27 Jun 2019 23:59:21 +0200 Subject: [PATCH 3/5] Using fast 2RDM s --- src/casscf/densities.irp.f | 146 ++++++--------------------------- src/casscf/test_two_rdm.irp.f | 10 +-- src/two_body_rdm/two_rdm.irp.f | 30 ++++--- 3 files changed, 46 insertions(+), 140 deletions(-) diff --git a/src/casscf/densities.irp.f b/src/casscf/densities.irp.f index 9b8dba78..7b243bb4 100644 --- a/src/casscf/densities.irp.f +++ b/src/casscf/densities.irp.f @@ -29,7 +29,9 @@ BEGIN_PROVIDER [real*8, P0tuvx, (n_act_orb,n_act_orb,n_act_orb,n_act_orb) ] ! END_DOC implicit none - integer :: t,u,v,x,mu,nu,istate,ispin,jspin,ihole,ipart,jhole,jpart + integer :: t,u,v,x + integer :: tt,uu,vv,xx + integer :: mu,nu,istate,ispin,jspin,ihole,ipart,jhole,jpart integer :: ierr real*8 :: phase1,phase11,phase12,phase2,phase21,phase22 integer :: nu1,nu2,nu11,nu12,nu21,nu22 @@ -43,125 +45,25 @@ BEGIN_PROVIDER [real*8, P0tuvx, (n_act_orb,n_act_orb,n_act_orb,n_act_orb) ] write(6,*) ' providing density matrix P0' endif - P0tuvx = 0.d0 - - ! first loop: we apply E_tu, once for D_tu, once for -P_tvvu - do mu=1,n_det - call det_extract(det_mu,mu,N_int) - do istate=1,n_states - cI_mu(istate)=psi_coef(mu,istate) - end do - do t=1,n_act_orb - ipart=list_act(t) - do u=1,n_act_orb - ihole=list_act(u) - ! apply E_tu - call det_copy(det_mu,det_mu_ex1,N_int) - call det_copy(det_mu,det_mu_ex2,N_int) - call do_spinfree_mono_excitation(det_mu,det_mu_ex1 & - ,det_mu_ex2,nu1,nu2,ihole,ipart,phase1,phase2,ierr1,ierr2) - ! det_mu_ex1 is in the list - if (nu1.ne.-1) then - do istate=1,n_states - term=cI_mu(istate)*psi_coef(nu1,istate)*phase1 - ! and we fill P0_tvvu - do v=1,n_act_orb - P0tuvx(t,v,v,u)-=term - end do - end do - end if - ! det_mu_ex2 is in the list - if (nu2.ne.-1) then - do istate=1,n_states - term=cI_mu(istate)*psi_coef(nu2,istate)*phase2 - do v=1,n_act_orb - P0tuvx(t,v,v,u)-=term - end do - end do - end if - end do - end do - end do - ! now we do the double excitation E_tu E_vx |0> - do mu=1,n_det - call det_extract(det_mu,mu,N_int) - do istate=1,n_states - cI_mu(istate)=psi_coef(mu,istate) - end do - do v=1,n_act_orb - ipart=list_act(v) - do x=1,n_act_orb - ihole=list_act(x) - ! apply E_vx - call det_copy(det_mu,det_mu_ex1,N_int) - call det_copy(det_mu,det_mu_ex2,N_int) - call do_spinfree_mono_excitation(det_mu,det_mu_ex1 & - ,det_mu_ex2,nu1,nu2,ihole,ipart,phase1,phase2,ierr1,ierr2) - ! we apply E_tu to the first resultant determinant, thus E_tu E_vx |0> - if (ierr1.eq.1) then - do t=1,n_act_orb - jpart=list_act(t) - do u=1,n_act_orb - jhole=list_act(u) - call det_copy(det_mu_ex1,det_mu_ex11,N_int) - call det_copy(det_mu_ex1,det_mu_ex12,N_int) - call do_spinfree_mono_excitation(det_mu_ex1,det_mu_ex11& - ,det_mu_ex12,nu11,nu12,jhole,jpart,phase11,phase12,ierr11,ierr12) - if (nu11.ne.-1) then - do istate=1,n_states - P0tuvx(t,u,v,x)+=cI_mu(istate)*psi_coef(nu11,istate)& - *phase11*phase1 - end do - end if - if (nu12.ne.-1) then - do istate=1,n_states - P0tuvx(t,u,v,x)+=cI_mu(istate)*psi_coef(nu12,istate)& - *phase12*phase1 - end do - end if - end do - end do - end if - - ! we apply E_tu to the second resultant determinant - if (ierr2.eq.1) then - do t=1,n_act_orb - jpart=list_act(t) - do u=1,n_act_orb - jhole=list_act(u) - call det_copy(det_mu_ex2,det_mu_ex21,N_int) - call det_copy(det_mu_ex2,det_mu_ex22,N_int) - call do_spinfree_mono_excitation(det_mu_ex2,det_mu_ex21& - ,det_mu_ex22,nu21,nu22,jhole,jpart,phase21,phase22,ierr21,ierr22) - if (nu21.ne.-1) then - do istate=1,n_states - P0tuvx(t,u,v,x)+=cI_mu(istate)*psi_coef(nu21,istate)& - *phase21*phase2 - end do - end if - if (nu22.ne.-1) then - do istate=1,n_states - P0tuvx(t,u,v,x)+=cI_mu(istate)*psi_coef(nu22,istate)& - *phase22*phase2 - end do - end if - end do - end do - end if - - end do - end do - end do - - ! we average by just dividing by the number of states - do x=1,n_act_orb - do v=1,n_act_orb - do u=1,n_act_orb - do t=1,n_act_orb - P0tuvx(t,u,v,x)*=0.5D0/dble(N_states) - end do - end do - end do - end do - + P0tuvx= 0.d0 + do istate=1,N_states + do x = 1, n_act_orb + xx = list_act(x) + do v = 1, n_act_orb + vv = list_act(v) + do u = 1, n_act_orb + uu = list_act(u) + do t = 1, n_act_orb + tt = list_act(t) + P0tuvx(t,u,v,x) = & + state_average_weight(istate) * & + ( two_rdm_alpha_beta_mo (tt,uu,vv,xx,istate) + & + two_rdm_alpha_alpha_mo(tt,uu,vv,xx,istate) + & + two_rdm_beta_beta_mo (tt,uu,vv,xx,istate) ) + enddo + enddo + enddo + enddo + enddo + END_PROVIDER diff --git a/src/casscf/test_two_rdm.irp.f b/src/casscf/test_two_rdm.irp.f index 562d15a6..9abe0aa0 100644 --- a/src/casscf/test_two_rdm.irp.f +++ b/src/casscf/test_two_rdm.irp.f @@ -8,11 +8,11 @@ program print_two_rdm double precision :: accu,twodm accu = 0.d0 - do i=1,mo_num - do j=1,mo_num - do k=1,mo_num - do l=1,mo_num - twodm = coussin_peter_two_rdm_mo(i,j,k,l,1) + do i=1,n_act_orb + do j=1,n_act_orb + do k=1,n_act_orb + do l=1,n_act_orb + twodm = coussin_peter_two_rdm_mo(list_act(i),list_act(j),list_act(k),list_act(l)) if(dabs(twodm - P0tuvx(i,j,k,l)).gt.thr)then print*,'' print*,'sum' diff --git a/src/two_body_rdm/two_rdm.irp.f b/src/two_body_rdm/two_rdm.irp.f index 1c299bba..a75a92cc 100644 --- a/src/two_body_rdm/two_rdm.irp.f +++ b/src/two_body_rdm/two_rdm.irp.f @@ -1,23 +1,27 @@ - - BEGIN_PROVIDER [double precision, coussin_peter_two_rdm_mo, (mo_num,mo_num,mo_num,mo_num,N_states)] +BEGIN_PROVIDER [double precision, coussin_peter_two_rdm_mo, (mo_num,mo_num,mo_num,mo_num)] implicit none BEGIN_DOC ! coussin_peter_two_rdm_mo(i,j,k,l) = the two rdm that peter wants for his CASSCF END_DOC - integer :: i,j,k,l - do l = 1, mo_num - do k = 1, mo_num - do j = 1, mo_num - do i = 1, mo_num - coussin_peter_two_rdm_mo(i,j,k,l,:) = 0.5d0 * (two_rdm_alpha_beta_mo(i,j,k,l,:) + two_rdm_alpha_beta_mo(i,j,k,l,:)) & - + two_rdm_alpha_alpha_mo(i,j,k,l,:) & - + two_rdm_beta_beta_mo(i,j,k,l,:) - enddo - enddo + integer :: i,j,k,l, istate + coussin_peter_two_rdm_mo = 0.d0 + do istate=1,N_states + do l = 1, mo_num + do k = 1, mo_num + do j = 1, mo_num + do i = 1, mo_num + coussin_peter_two_rdm_mo(i,j,k,l) = & + state_average_weight(istate) * & + ( two_rdm_alpha_beta_mo(i,j,k,l,istate) + & + two_rdm_alpha_alpha_mo(i,j,k,l,istate)+ & + two_rdm_beta_beta_mo(i,j,k,l,istate) ) + enddo + enddo + enddo enddo enddo - END_PROVIDER +END_PROVIDER BEGIN_PROVIDER [double precision, two_rdm_alpha_beta_mo, (mo_num,mo_num,mo_num,mo_num,N_states)] From a4d2e39978ecb4802fee53b80e9d525e8e010900 Mon Sep 17 00:00:00 2001 From: Anthony Scemama Date: Fri, 28 Jun 2019 00:04:12 +0200 Subject: [PATCH 4/5] Minor fix --- src/cipsi/cipsi.irp.f | 1 + src/cipsi/stochastic_cipsi.irp.f | 1 + 2 files changed, 2 insertions(+) diff --git a/src/cipsi/cipsi.irp.f b/src/cipsi/cipsi.irp.f index 7e292d6e..ba922c49 100644 --- a/src/cipsi/cipsi.irp.f +++ b/src/cipsi/cipsi.irp.f @@ -13,6 +13,7 @@ subroutine run_cipsi rss = memory_of_double(N_states)*4.d0 call check_mem(rss,irp_here) + N_iter = 1 allocate (pt2(N_states), zeros(N_states), rpt2(N_states), norm(N_states), variance(N_states)) double precision :: hf_energy_ref diff --git a/src/cipsi/stochastic_cipsi.irp.f b/src/cipsi/stochastic_cipsi.irp.f index ae2b7519..4f968ef7 100644 --- a/src/cipsi/stochastic_cipsi.irp.f +++ b/src/cipsi/stochastic_cipsi.irp.f @@ -12,6 +12,7 @@ subroutine run_stochastic_cipsi double precision, external :: memory_of_double PROVIDE H_apply_buffer_allocated N_generators_bitmask + N_iter = 1 threshold_generators = 1.d0 SOFT_TOUCH threshold_generators From d742bdd655a28648652dc9ba45ea96248a32ce11 Mon Sep 17 00:00:00 2001 From: Anthony Scemama Date: Fri, 28 Jun 2019 00:06:51 +0200 Subject: [PATCH 5/5] Cleaning --- src/casscf/test_two_rdm.irp.f | 30 ------------------------------ src/two_body_rdm/two_rdm.irp.f | 26 -------------------------- 2 files changed, 56 deletions(-) delete mode 100644 src/casscf/test_two_rdm.irp.f diff --git a/src/casscf/test_two_rdm.irp.f b/src/casscf/test_two_rdm.irp.f deleted file mode 100644 index 9abe0aa0..00000000 --- a/src/casscf/test_two_rdm.irp.f +++ /dev/null @@ -1,30 +0,0 @@ -program print_two_rdm - implicit none - integer :: i,j,k,l - read_wf = .True. - TOUCH read_wf - - double precision, parameter :: thr = 1.d-15 - - double precision :: accu,twodm - accu = 0.d0 - do i=1,n_act_orb - do j=1,n_act_orb - do k=1,n_act_orb - do l=1,n_act_orb - twodm = coussin_peter_two_rdm_mo(list_act(i),list_act(j),list_act(k),list_act(l)) - if(dabs(twodm - P0tuvx(i,j,k,l)).gt.thr)then - print*,'' - print*,'sum' - write(*,'(3X,4(I2,X),3(F16.13,X))'), i, j, k, l, twodm,P0tuvx(i,j,k,l),dabs(twodm - P0tuvx(i,j,k,l)) - print*,'' - endif - accu += dabs(twodm - P0tuvx(i,j,k,l)) - enddo - enddo - enddo - enddo - print*,'accu = ',accu - print*,' ',accu / dble(mo_num**4) - -end diff --git a/src/two_body_rdm/two_rdm.irp.f b/src/two_body_rdm/two_rdm.irp.f index a75a92cc..89eecdcc 100644 --- a/src/two_body_rdm/two_rdm.irp.f +++ b/src/two_body_rdm/two_rdm.irp.f @@ -1,29 +1,3 @@ -BEGIN_PROVIDER [double precision, coussin_peter_two_rdm_mo, (mo_num,mo_num,mo_num,mo_num)] - implicit none - BEGIN_DOC - ! coussin_peter_two_rdm_mo(i,j,k,l) = the two rdm that peter wants for his CASSCF - END_DOC - integer :: i,j,k,l, istate - coussin_peter_two_rdm_mo = 0.d0 - do istate=1,N_states - do l = 1, mo_num - do k = 1, mo_num - do j = 1, mo_num - do i = 1, mo_num - coussin_peter_two_rdm_mo(i,j,k,l) = & - state_average_weight(istate) * & - ( two_rdm_alpha_beta_mo(i,j,k,l,istate) + & - two_rdm_alpha_alpha_mo(i,j,k,l,istate)+ & - two_rdm_beta_beta_mo(i,j,k,l,istate) ) - enddo - enddo - enddo - enddo - enddo - -END_PROVIDER - - BEGIN_PROVIDER [double precision, two_rdm_alpha_beta_mo, (mo_num,mo_num,mo_num,mo_num,N_states)] &BEGIN_PROVIDER [double precision, two_rdm_alpha_alpha_mo, (mo_num,mo_num,mo_num,mo_num,N_states)] &BEGIN_PROVIDER [double precision, two_rdm_beta_beta_mo, (mo_num,mo_num,mo_num,mo_num,N_states)]