diff --git a/src/casscf_cipsi/gradient.irp.f b/src/casscf_cipsi/gradient.irp.f index 961d260d..cd608fb3 100644 --- a/src/casscf_cipsi/gradient.irp.f +++ b/src/casscf_cipsi/gradient.irp.f @@ -14,8 +14,8 @@ END_PROVIDER implicit none n_c_a_prov = n_core_inact_orb * n_act_orb n_c_v_prov = n_core_inact_orb * n_virt_orb - n_a_v_prov = n_act_orb * n_virt_orb - END_PROVIDER + n_a_v_prov = n_act_orb * n_virt_orb + END_PROVIDER BEGIN_PROVIDER [integer, excit, (2,nMonoEx)] &BEGIN_PROVIDER [character*3, excit_class, (nMonoEx)] @@ -28,7 +28,7 @@ END_PROVIDER BEGIN_DOC ! a list of the orbitals involved in the excitation END_DOC - + implicit none integer :: i,t,a,ii,tt,aa,indx,indx_tmp indx=0 @@ -48,7 +48,7 @@ END_PROVIDER mat_idx_c_a(ii,tt) = indx end do end do - + indx_tmp = 0 do ii=1,n_core_inact_orb i=list_core_inact(ii) @@ -61,11 +61,11 @@ END_PROVIDER indx_tmp += 1 list_idx_c_v(1,indx_tmp) = indx list_idx_c_v(2,indx_tmp) = ii - list_idx_c_v(3,indx_tmp) = aa + list_idx_c_v(3,indx_tmp) = aa mat_idx_c_v(ii,aa) = indx end do end do - + indx_tmp = 0 do tt=1,n_act_orb t=list_act(tt) @@ -82,7 +82,7 @@ END_PROVIDER mat_idx_a_v(tt,aa) = indx end do end do - + if (bavard) then write(6,*) ' Filled the table of the Monoexcitations ' do indx=1,nMonoEx @@ -90,7 +90,7 @@ END_PROVIDER ,excit(2,indx),' ',excit_class(indx) end do end if - + END_PROVIDER BEGIN_PROVIDER [real*8, gradvec2, (nMonoEx)] @@ -104,7 +104,7 @@ END_PROVIDER implicit none integer :: i,t,a,indx real*8 :: gradvec_it,gradvec_ia,gradvec_ta - + indx=0 norm_grad_vec2_tab = 0.d0 do i=1,n_core_inact_orb @@ -114,7 +114,7 @@ END_PROVIDER norm_grad_vec2_tab(1) += gradvec2(indx)*gradvec2(indx) end do end do - + do i=1,n_core_inact_orb do a=1,n_virt_orb indx+=1 @@ -122,7 +122,7 @@ END_PROVIDER norm_grad_vec2_tab(2) += gradvec2(indx)*gradvec2(indx) end do end do - + do t=1,n_act_orb do a=1,n_virt_orb indx+=1 @@ -130,7 +130,7 @@ END_PROVIDER norm_grad_vec2_tab(3) += gradvec2(indx)*gradvec2(indx) end do end do - + norm_grad_vec2=0.d0 do indx=1,nMonoEx norm_grad_vec2+=gradvec2(indx)*gradvec2(indx) @@ -144,7 +144,7 @@ END_PROVIDER write(6,*) ' Norm of the orbital gradient (via D, P and integrals): ', norm_grad_vec2 write(6,*) endif - + END_PROVIDER real*8 function gradvec_it(i,t) @@ -154,23 +154,30 @@ real*8 function gradvec_it(i,t) END_DOC implicit none integer :: i,t - + integer :: ii,tt,v,vv,x,y integer :: x3,y3 double precision :: bielec_PQxx_no - + ii=list_core_inact(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 ! active - vv=list_act(v) - do x=1,n_act_orb ! active - x3=x+n_core_inact_orb ! list_act(x) - do y=1,n_act_orb ! active - y3=y+n_core_inact_orb ! list_act(y) + do y=1,n_act_orb ! active +! y3=y+n_core_inact_orb ! list_act(y) + do x=1,n_act_orb ! active +! x3=x+n_core_inact_orb ! list_act(x) + do v=1,n_act_orb ! active + vv=list_act(v) ! Gamma(2) a a a a 1/r12 i a a a - gradvec_it-=2.D0*P0tuvx_no(t,v,x,y)*bielec_PQxx_no(ii,vv,x3,y3) +! gradvec_it-=2.D0*P0tuvx_no(t,v,x,y)*bielec_PQxx_no(ii,vv,x3,y3) + integer :: ichol + double precision :: tmp + tmp = 0.d0 + do ichol=1,cholesky_mo_num + tmp = tmp + cholesky_no_total_transp(ichol,vv,ii) * cholesky_no_total_transp(ichol,list_act(x),list_act(y)) + enddo + gradvec_it = gradvec_it - 2.D0*P0tuvx_no(t,v,x,y)*tmp end do end do end do @@ -183,12 +190,12 @@ real*8 function gradvec_ia(i,a) END_DOC implicit none integer :: i,a,ii,aa - + ii=list_core_inact(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) @@ -198,7 +205,7 @@ real*8 function gradvec_ta(t,a) 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 @@ -211,6 +218,6 @@ real*8 function gradvec_ta(t,a) end do end do gradvec_ta*=2.D0 - + end function gradvec_ta diff --git a/src/casscf_cipsi/hessian.irp.f b/src/casscf_cipsi/hessian.irp.f index 1ee073d2..1bcb64b0 100644 --- a/src/casscf_cipsi/hessian.irp.f +++ b/src/casscf_cipsi/hessian.irp.f @@ -11,13 +11,14 @@ 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 double precision :: bielec_pqxx_no,bielec_pxxq_no - + ii=list_core_inact(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)) & + 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_no(tt,i,i,tt)-bielec_pqxx_no(tt,tt,i,i)) term-=2.D0*occnum(tt)*(3.D0*bielec_pxxq_no(tt,i,i,tt) & @@ -83,10 +84,10 @@ real*8 function hessmat_itju(i,t,j,u) 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) @@ -97,7 +98,7 @@ 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 double precision :: bielec_pqxx_no,bielec_pxxq_no - + ! it/ja ii=list_core_inact(i) tt=list_act(t) @@ -120,7 +121,7 @@ real*8 function hessmat_itja(i,t,j,a) end if term*=2.D0 hessmat_itja=term - + end function hessmat_itja real*8 function hessmat_itua(i,t,u,a) @@ -131,7 +132,7 @@ 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 double precision :: bielec_pqxx_no,bielec_pxxq_no - + ii=list_core_inact(i) tt=list_act(t) t3=t+n_core_inact_orb @@ -162,7 +163,7 @@ real*8 function hessmat_itua(i,t,u,a) end if term*=2.D0 hessmat_itua=term - + end function hessmat_itua real*8 function hessmat_iajb(i,a,j,b) @@ -173,7 +174,7 @@ real*8 function hessmat_iajb(i,a,j,b) integer :: i,a,j,b,ii,aa,jj,bb real*8 :: term double precision :: bielec_pqxx_no,bielec_pxxq_no - + ii=list_core_inact(i) aa=list_virt(a) if (i.eq.j) then @@ -199,7 +200,7 @@ real*8 function hessmat_iajb(i,a,j,b) end if term*=2.D0 hessmat_iajb=term - + end function hessmat_iajb real*8 function hessmat_iatb(i,a,t,b) @@ -210,7 +211,7 @@ 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 double precision :: bielec_pqxx_no,bielec_pxxq_no - + ii=list_core_inact(i) aa=list_virt(a) tt=list_act(t) @@ -231,7 +232,7 @@ real*8 function hessmat_iatb(i,a,t,b) end if term*=2.D0 hessmat_iatb=term - + end function hessmat_iatb real*8 function hessmat_taub(t,a,u,b) @@ -240,83 +241,186 @@ real*8 function hessmat_taub(t,a,u,b) 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 + integer :: v3,x3, ichol + real*8 :: term,t1,t2,t3, tmp double precision :: bielec_pqxx_no,bielec_pxxq_no - + + double precision, allocatable :: tmp1(:), tmp2(:,:) + allocate(tmp1(n_act_orb)) + allocate(tmp2(n_act_orb,n_act_orb)) + tt=list_act(t) aa=list_virt(a) + if (t == u) then if (a == b) then ! ta/ta - t1=occnum(tt)*Fipq(aa,aa) + t1=occnum(tt)*Fipq(aa,aa) - occnum(tt)*Fipq(tt,tt) + t2=0.D0 - t3=0.D0 - t1-=occnum(tt)*Fipq(tt,tt) +! do x=1,n_act_orb +! x3=x+n_core_inact_orb +! do v=1,n_act_orb +! v3=v+n_core_inact_orb +! tmp = 0.d0 +! do ichol = 1, cholesky_mo_num +! tmp = tmp + cholesky_no_total_transp(ichol,aa,aa) * cholesky_no_total_transp(ichol,v3,x3) +! enddo +! t2 = t2 + 2.D0*P0tuvx_no(t,t,v,x)*tmp +! enddo +! enddo + + do x=1,n_act_orb + x3=x+n_core_inact_orb + call dgemv('T', cholesky_mo_num, n_act_orb, 2.d0, & + cholesky_no_total_transp(1,n_core_inact_orb+1,x3), cholesky_mo_num, & + cholesky_no_total_transp(1,aa,aa), 1, 0.d0, & + tmp1, 1) + do v=1,n_act_orb + t2 = t2 + P0tuvx_no(t,t,v,x)*tmp1(v) + enddo + enddo +! do v=1,n_act_orb +! v3=v+n_core_inact_orb +! do x=1,n_act_orb +! x3=x+n_core_inact_orb +! tmp = 0.d0 +! do ichol = 1, cholesky_mo_num +! tmp = tmp + cholesky_no_total_transp(ichol,aa,x3) * cholesky_no_total_transp(ichol,v3,aa) +! enddo +! t2 = t2 + 2.d0*(P0tuvx_no(t,x,v,t)+P0tuvx_no(t,x,t,v))*tmp +! end do +! end do + call dgemm('T','N', n_act_orb, n_act_orb, cholesky_mo_num, 2.d0, & + cholesky_no_total_transp(1,n_core_inact_orb+1,aa), cholesky_mo_num, & + cholesky_no_total_transp(1,n_core_inact_orb+1,aa), cholesky_mo_num, 0.d0, & + tmp2, n_act_orb) do v=1,n_act_orb - vv=list_act(v) - v3=v+n_core_inact_orb do x=1,n_act_orb - xx=list_act(x) - x3=x+n_core_inact_orb - t2+=2.D0*(P0tuvx_no(t,t,v,x)*bielec_pqxx_no(aa,aa,v3,x3) & - +(P0tuvx_no(t,x,v,t)+P0tuvx_no(t,x,t,v))* & - bielec_pxxq_no(aa,x3,v3,aa)) - do y=1,n_act_orb - t3-=2.D0*P0tuvx_no(t,v,x,y)*bielecCI_no(t,v,y,xx) + t2 = t2 + P0tuvx_no(t,x,v,t)*tmp2(x,v) + P0tuvx_no(t,x,t,v)*tmp2(x,v) + enddo + enddo + + t3=0.D0 + do x=1,n_act_orb + xx=list_act(x) + do y=1,n_act_orb + do v=1,n_act_orb + t3 = t3 - P0tuvx_no(t,v,x,y)*bielecCI_no(v,t,y,xx) end do end do end do - term=t1+t2+t3 + term=t1+t2+t3*2.d0 + 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_inact_orb +! do x=1,n_act_orb +! xx=list_act(x) +! x3=x+n_core_inact_orb +! term+=2.D0*P0tuvx_no(t,t,v,x)*bielec_pqxx_no(aa,bb,v3,x3) +! end do +! end do + do x=1,n_act_orb + x3=x+n_core_inact_orb + call dgemv('T', cholesky_mo_num, n_act_orb, 2.d0, & + cholesky_no_total_transp(1,n_core_inact_orb+1,x3), cholesky_mo_num, & + cholesky_no_total_transp(1,aa,bb), 1, 0.d0, & + tmp1, 1) + do v=1,n_act_orb + term = term + P0tuvx_no(t,t,v,x)*tmp1(v) + enddo + enddo + +! do v=1,n_act_orb +! vv=list_act(v) +! v3=v+n_core_inact_orb +! do x=1,n_act_orb +! xx=list_act(x) +! x3=x+n_core_inact_orb +! term+=2.d0*(P0tuvx_no(t,x,v,t)+P0tuvx_no(t,x,t,v))*bielec_pxxq_no(aa,x3,v3,bb) +! end do +! end do + call dgemm('T','N', n_act_orb, n_act_orb, cholesky_mo_num, 2.d0, & + cholesky_no_total_transp(1,n_core_inact_orb+1,aa), cholesky_mo_num, & + cholesky_no_total_transp(1,n_core_inact_orb+1,bb), cholesky_mo_num, 0.d0, & + tmp2, n_act_orb) do v=1,n_act_orb - vv=list_act(v) - v3=v+n_core_inact_orb do x=1,n_act_orb - xx=list_act(x) - x3=x+n_core_inact_orb - term+=2.D0*(P0tuvx_no(t,t,v,x)*bielec_pqxx_no(aa,bb,v3,x3) & - +(P0tuvx_no(t,x,v,t)+P0tuvx_no(t,x,t,v)) & - *bielec_pxxq_no(aa,x3,v3,bb)) - end do - end do + term = term + P0tuvx_no(t,x,v,t)*tmp2(x,v) + P0tuvx_no(t,x,t,v)*tmp2(x,v) + enddo + enddo + 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_inact_orb +! do x=1,n_act_orb +! xx=list_act(x) +! x3=x+n_core_inact_orb +! term+=2.D0*(P0tuvx_no(t,u,v,x)*bielec_pqxx_no(aa,bb,v3,x3) +! end do +! end do + do x=1,n_act_orb + x3=x+n_core_inact_orb + call dgemv('T', cholesky_mo_num, n_act_orb, 2.d0, & + cholesky_no_total_transp(1,n_core_inact_orb+1,x3), cholesky_mo_num, & + cholesky_no_total_transp(1,aa,bb), 1, 0.d0, & + tmp1, 1) + do v=1,n_act_orb + term = term + P0tuvx_no(t,u,v,x)*tmp1(v) + enddo + enddo + +! do v=1,n_act_orb +! vv=list_act(v) +! v3=v+n_core_inact_orb +! do x=1,n_act_orb +! xx=list_act(x) +! x3=x+n_core_inact_orb +! term+=2.D0*(P0tuvx_no(t,x,v,u)+P0tuvx_no(t,x,u,v))*bielec_pxxq_no(aa,x3,v3,bb) +! end do +! end do + call dgemm('T','N', n_act_orb, n_act_orb, cholesky_mo_num, 2.d0, & + cholesky_no_total_transp(1,n_core_inact_orb+1,aa), cholesky_mo_num, & + cholesky_no_total_transp(1,n_core_inact_orb+1,bb), cholesky_mo_num, 0.d0, & + tmp2, n_act_orb) do v=1,n_act_orb - vv=list_act(v) - v3=v+n_core_inact_orb do x=1,n_act_orb - xx=list_act(x) - x3=x+n_core_inact_orb - term+=2.D0*(P0tuvx_no(t,u,v,x)*bielec_pqxx_no(aa,bb,v3,x3) & - +(P0tuvx_no(t,x,v,u)+P0tuvx_no(t,x,u,v)) & - *bielec_pxxq_no(aa,x3,v3,bb)) - end do - end do + term = term + P0tuvx_no(t,x,v,u)*tmp2(x,v)+P0tuvx_no(t,x,u,v)*tmp2(x,v) + enddo + enddo + 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 y=1,n_act_orb do x=1,n_act_orb - term-=P0tuvx_no(t,v,x,y)*bielecCI_no(x,y,v,uu) - term-=P0tuvx_no(u,v,x,y)*bielecCI_no(x,y,v,tt) + term = term - P0tuvx_no(t,v,x,y)*bielecCI_no(x,y,v,uu) & + - P0tuvx_no(u,v,x,y)*bielecCI_no(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)] @@ -326,7 +430,7 @@ BEGIN_PROVIDER [real*8, hessdiag, (nMonoEx)] implicit none integer :: i,t,a,indx,indx_shift real*8 :: hessmat_itju,hessmat_iajb,hessmat_taub - + !$OMP PARALLEL DEFAULT(NONE) & !$OMP SHARED(hessdiag,n_core_inact_orb,n_act_orb,n_virt_orb,nMonoEx) & !$OMP PRIVATE(i,indx,t,a,indx_shift) @@ -339,9 +443,9 @@ BEGIN_PROVIDER [real*8, hessdiag, (nMonoEx)] end do end do !$OMP END DO NOWAIT - + indx_shift = n_core_inact_orb*n_act_orb - !$OMP DO + !$OMP DO do a=1,n_virt_orb do i=1,n_core_inact_orb indx = a + (i-1)*n_virt_orb + indx_shift @@ -349,9 +453,9 @@ BEGIN_PROVIDER [real*8, hessdiag, (nMonoEx)] end do end do !$OMP END DO NOWAIT - + indx_shift += n_core_inact_orb*n_virt_orb - !$OMP DO + !$OMP DO do a=1,n_virt_orb do t=1,n_act_orb indx = a + (t-1)*n_virt_orb + indx_shift @@ -360,7 +464,7 @@ BEGIN_PROVIDER [real*8, hessdiag, (nMonoEx)] end do !$OMP END DO !$OMP END PARALLEL - + END_PROVIDER @@ -377,7 +481,7 @@ BEGIN_PROVIDER [double precision, hessmat, (nMonoEx,nMonoEx)] real*8 :: hessmat_taub ! c-a c-v a-v ! c-a | X X X - ! c-v | X X + ! c-v | X X ! a-v | X provide all_mo_integrals @@ -390,12 +494,12 @@ BEGIN_PROVIDER [double precision, hessmat, (nMonoEx,nMonoEx)] !$OMP DO !!!! < Core-active| H |Core-active > - ! Core-active excitations + ! Core-active excitations do indx_tmp = 1, n_c_a_prov indx = list_idx_c_a(1,indx_tmp) i = list_idx_c_a(2,indx_tmp) t = list_idx_c_a(3,indx_tmp) - ! Core-active excitations + ! Core-active excitations do j = 1, n_core_inact_orb if (i.eq.j) then ustart=t @@ -418,12 +522,12 @@ BEGIN_PROVIDER [double precision, hessmat, (nMonoEx,nMonoEx)] !$OMP DO !!!! < Core-active| H |Core-VIRTUAL > - ! Core-active excitations + ! Core-active excitations do indx_tmp = 1, n_c_a_prov indx = list_idx_c_a(1,indx_tmp) i = list_idx_c_a(2,indx_tmp) t = list_idx_c_a(3,indx_tmp) - ! Core-VIRTUAL excitations + ! Core-VIRTUAL excitations do jndx_tmp = 1, n_c_v_prov jndx = list_idx_c_v(1,jndx_tmp) j = list_idx_c_v(2,jndx_tmp) @@ -441,12 +545,12 @@ BEGIN_PROVIDER [double precision, hessmat, (nMonoEx,nMonoEx)] !$OMP DO !!!! < Core-active| H |ACTIVE-VIRTUAL > - ! Core-active excitations + ! Core-active excitations do indx_tmp = 1, n_c_a_prov indx = list_idx_c_a(1,indx_tmp) i = list_idx_c_a(2,indx_tmp) t = list_idx_c_a(3,indx_tmp) - ! ACTIVE-VIRTUAL excitations + ! ACTIVE-VIRTUAL excitations do jndx_tmp = 1, n_a_v_prov jndx = list_idx_a_v(1,jndx_tmp) u = list_idx_a_v(2,jndx_tmp) @@ -466,12 +570,12 @@ BEGIN_PROVIDER [double precision, hessmat, (nMonoEx,nMonoEx)] !$OMP PRIVATE(indx_tmp,indx,i,a,j,b,bstart,jndx) !$OMP DO !!!!! < Core-VIRTUAL | H |Core-VIRTUAL > - ! Core-VIRTUAL excitations + ! Core-VIRTUAL excitations do indx_tmp = 1, n_c_v_prov indx = list_idx_c_v(1,indx_tmp) i = list_idx_c_v(2,indx_tmp) a = list_idx_c_v(3,indx_tmp) - ! Core-VIRTUAL excitations + ! Core-VIRTUAL excitations do j = 1, n_core_inact_orb if (i.eq.j) then bstart=a @@ -485,7 +589,7 @@ BEGIN_PROVIDER [double precision, hessmat, (nMonoEx,nMonoEx)] enddo enddo enddo - + !$OMP END DO NOWAIT !$OMP END PARALLEL endif @@ -496,12 +600,12 @@ BEGIN_PROVIDER [double precision, hessmat, (nMonoEx,nMonoEx)] !$OMP DO !!!! < Core-VIRTUAL | H |Active-VIRTUAL > - ! Core-VIRTUAL excitations + ! Core-VIRTUAL excitations do indx_tmp = 1, n_c_v_prov indx = list_idx_c_v(1,indx_tmp) i = list_idx_c_v(2,indx_tmp) a = list_idx_c_v(3,indx_tmp) - ! Active-VIRTUAL excitations + ! Active-VIRTUAL excitations do jndx_tmp = 1, n_a_v_prov jndx = list_idx_a_v(1,jndx_tmp) t = list_idx_a_v(2,jndx_tmp) @@ -520,12 +624,12 @@ BEGIN_PROVIDER [double precision, hessmat, (nMonoEx,nMonoEx)] !$OMP DO !!!! < Active-VIRTUAL | H |Active-VIRTUAL > - ! Active-VIRTUAL excitations + ! Active-VIRTUAL excitations do indx_tmp = 1, n_a_v_prov indx = list_idx_a_v(1,indx_tmp) t = list_idx_a_v(2,indx_tmp) a = list_idx_a_v(3,indx_tmp) - ! Active-VIRTUAL excitations + ! Active-VIRTUAL excitations do u=t,n_act_orb if (t.eq.u) then bstart=a @@ -542,4 +646,4 @@ BEGIN_PROVIDER [double precision, hessmat, (nMonoEx,nMonoEx)] !$OMP END DO NOWAIT !$OMP END PARALLEL -END_PROVIDER +END_PROVIDER diff --git a/src/casscf_cipsi/mcscf_fock.irp.f b/src/casscf_cipsi/mcscf_fock.irp.f index 82b710a7..87cc062c 100644 --- a/src/casscf_cipsi/mcscf_fock.irp.f +++ b/src/casscf_cipsi/mcscf_fock.irp.f @@ -3,37 +3,53 @@ BEGIN_PROVIDER [real*8, Fipq, (mo_num,mo_num) ] ! the inactive Fock matrix, in molecular orbitals END_DOC implicit none - integer :: p,q,k,kk,t,tt,u,uu - double precision :: bielec_pxxq_no, bielec_pqxx_no - + integer :: i,p,q,k,kk,t,tt,u,uu + do q=1,mo_num do p=1,mo_num Fipq(p,q)=one_ints_no(p,q) end do end do - + ! the inactive Fock matrix do k=1,n_core_inact_orb kk=list_core_inact(k) - do q=1,mo_num - do p=1,mo_num - Fipq(p,q)+=2.D0*bielec_pqxx_no(p,q,k,k) -bielec_pxxq_no(p,k,k,q) - end do - end do +! do q=1,mo_num +! do p=1,mo_num +! do i=1,cholesky_mo_num +! Fipq(p,q) = Fipq(p,q) + 2.d0* cholesky_no_total_transp(i,p,q) * cholesky_no_total_transp(i,kk,kk) +! enddo +! end do +! end do + call dgemv('T', cholesky_mo_num, mo_num*mo_num, 2.d0, & + cholesky_no_total_transp, cholesky_mo_num, & + cholesky_no_total_transp(1,kk,kk), 1, 1.d0, & + Fipq, 1) + +! do q=1,mo_num +! do p=1,mo_num +! do i=1,cholesky_mo_num +! Fipq(p,q) = Fipq(p,q) - cholesky_no_total_transp(i,p,kk) * cholesky_no_total_transp(i,kk,q) +! enddo +! end do +! end do + call dgemm('T','N', mo_num, mo_num, cholesky_mo_num, -1.d0, & + cholesky_no_total_transp(1,1,kk), cholesky_mo_num, & + cholesky_no_total_transp(1,kk,1), cholesky_mo_num*mo_num, 1.d0, & + Fipq, mo_num) end do - + if (bavard) then - integer :: i 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,*) end if - - + + END_PROVIDER - - + + BEGIN_PROVIDER [real*8, Fapq, (mo_num,mo_num) ] BEGIN_DOC ! the active active Fock matrix, in molecular orbitals @@ -45,27 +61,42 @@ BEGIN_PROVIDER [real*8, Fapq, (mo_num,mo_num) ] END_DOC implicit none integer :: p,q,k,kk,t,tt,u,uu - double precision :: bielec_pxxq_no, bielec_pqxx_no - + Fapq = 0.d0 - + ! the active Fock matrix, D0tu is diagonal do t=1,n_act_orb tt=list_act(t) - do q=1,mo_num - do p=1,mo_num - Fapq(p,q)+=occnum(tt) & - *(bielec_pqxx_no(p,q,tt,tt)-0.5D0*bielec_pxxq_no(p,tt,tt,q)) - end do - end do +! do q=1,mo_num +! do p=1,mo_num +! do i=1,cholesky_mo_num +! Fapq(p,q) = Fapq(p,q) + occnum(tt)* cholesky_no_total_transp(i,p,q) * cholesky_no_total_transp(i,tt,tt) +! enddo +! end do +! end do + call dgemv('T', cholesky_mo_num, mo_num*mo_num, occnum(tt), & + cholesky_no_total_transp, cholesky_mo_num, & + cholesky_no_total_transp(1,tt,tt), 1, 1.d0, & + Fapq, 1) +! do q=1,mo_num +! do p=1,mo_num +! do i=1,cholesky_mo_num +! Fapq(p,q) = Fapq(p,q) - 0.5d0*occnum(tt)*cholesky_no_total_transp(i,p,tt) * cholesky_no_total_transp(i,tt,q) +! enddo +! end do +! end do + call dgemm('T','N', mo_num, mo_num, cholesky_mo_num, -0.5d0*occnum(tt), & + cholesky_no_total_transp(1,1,tt), cholesky_mo_num, & + cholesky_no_total_transp(1,tt,1), cholesky_mo_num*mo_num, 1.d0, & + Fapq, mo_num) 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) @@ -75,35 +106,35 @@ BEGIN_PROVIDER [real*8, Fapq, (mo_num,mo_num) ] write(6,'(5(i3,F12.5))') (i,Fapq(i,i),i=1,mo_num) write(6,*) end if - - + + END_PROVIDER - - BEGIN_PROVIDER [ double precision, mcscf_fock_alpha_ao, (ao_num, ao_num)] -&BEGIN_PROVIDER [ double precision, mcscf_fock_beta_ao, (ao_num, ao_num)] + + BEGIN_PROVIDER [ double precision, mcscf_fock_alpha_ao, (ao_num, ao_num)] +&BEGIN_PROVIDER [ double precision, mcscf_fock_beta_ao, (ao_num, ao_num)] implicit none BEGIN_DOC - ! mcscf_fock_alpha_ao are set to usual Fock like operator but computed with the MCSCF densities on the AO basis + ! mcscf_fock_alpha_ao are set to usual Fock like operator but computed with the MCSCF densities on the AO basis END_DOC SCF_density_matrix_ao_alpha = D0tu_alpha_ao SCF_density_matrix_ao_beta = D0tu_beta_ao - soft_touch SCF_density_matrix_ao_alpha SCF_density_matrix_ao_beta + soft_touch SCF_density_matrix_ao_alpha SCF_density_matrix_ao_beta mcscf_fock_beta_ao = fock_matrix_ao_beta mcscf_fock_alpha_ao = fock_matrix_ao_alpha -END_PROVIDER +END_PROVIDER - BEGIN_PROVIDER [ double precision, mcscf_fock_alpha_mo, (mo_num, mo_num)] -&BEGIN_PROVIDER [ double precision, mcscf_fock_beta_mo, (mo_num, mo_num)] + BEGIN_PROVIDER [ double precision, mcscf_fock_alpha_mo, (mo_num, mo_num)] +&BEGIN_PROVIDER [ double precision, mcscf_fock_beta_mo, (mo_num, mo_num)] implicit none BEGIN_DOC - ! Mo_mcscf_fock_alpha are set to usual Fock like operator but computed with the MCSCF densities on the MO basis + ! Mo_mcscf_fock_alpha are set to usual Fock like operator but computed with the MCSCF densities on the MO basis END_DOC call ao_to_mo(mcscf_fock_alpha_ao,ao_num,mcscf_fock_alpha_mo,mo_num) call ao_to_mo(mcscf_fock_beta_ao,ao_num,mcscf_fock_beta_mo,mo_num) -END_PROVIDER +END_PROVIDER BEGIN_PROVIDER [ double precision, mcscf_fock_mo, (mo_num,mo_num) ] &BEGIN_PROVIDER [ double precision, mcscf_fock_diag_mo, (mo_num)] @@ -118,13 +149,13 @@ END_PROVIDER ! |-----------------------| ! | Fcv | F^a | Rvv | ! - ! C: Core, O: Open, V: Virtual - ! + ! C: Core, O: Open, V: Virtual + ! ! Rcc = Acc Fcc^a + Bcc Fcc^b ! Roo = Aoo Foo^a + Boo Foo^b ! Rvv = Avv Fvv^a + Bvv Fvv^b ! Fcv = (F^a + F^b)/2 - ! + ! ! F^a: Fock matrix alpha (MO), F^b: Fock matrix beta (MO) ! A,B: Coupling parameters ! @@ -133,7 +164,7 @@ END_PROVIDER ! cc oo vv ! A -0.5 0.5 1.5 ! B 1.5 0.5 -0.5 - ! + ! END_DOC integer :: i,j,n if (elec_alpha_num == elec_beta_num) then @@ -194,4 +225,4 @@ END_PROVIDER do i = 1, mo_num mcscf_fock_diag_mo(i) = mcscf_fock_mo(i,i) enddo -END_PROVIDER +END_PROVIDER diff --git a/src/casscf_cipsi/natorb.irp.f b/src/casscf_cipsi/natorb.irp.f index 9ce90304..6376308d 100644 --- a/src/casscf_cipsi/natorb.irp.f +++ b/src/casscf_cipsi/natorb.irp.f @@ -72,84 +72,27 @@ BEGIN_PROVIDER [real*8, P0tuvx_no, (n_act_orb,n_act_orb,n_act_orb,n_act_orb)] BEGIN_DOC ! 4-index transformation of 2part matrices END_DOC - integer :: i,j,k,l,p,q - real*8 :: d(n_act_orb) - ! index per index - ! first quarter - P0tuvx_no(:,:,:,:) = P0tuvx(:,:,:,:) + double precision, allocatable :: tmp(:,:,:,:) + allocate(tmp(n_act_orb,n_act_orb,n_act_orb,n_act_orb)) - do j=1,n_act_orb - do k=1,n_act_orb - do l=1,n_act_orb - do p=1,n_act_orb - d(p)=0.D0 - end do - do p=1,n_act_orb - do q=1,n_act_orb - d(p)+=P0tuvx_no(q,j,k,l)*natorbsCI(q,p) - end do - end do - do p=1,n_act_orb - P0tuvx_no(p,j,k,l)=d(p) - end do - end do - end do - end do - ! 2nd quarter - do j=1,n_act_orb - do k=1,n_act_orb - do l=1,n_act_orb - do p=1,n_act_orb - d(p)=0.D0 - end do - do p=1,n_act_orb - do q=1,n_act_orb - d(p)+=P0tuvx_no(j,q,k,l)*natorbsCI(q,p) - end do - end do - do p=1,n_act_orb - P0tuvx_no(j,p,k,l)=d(p) - end do - end do - end do - end do - ! 3rd quarter - do j=1,n_act_orb - do k=1,n_act_orb - do l=1,n_act_orb - do p=1,n_act_orb - d(p)=0.D0 - end do - do p=1,n_act_orb - do q=1,n_act_orb - d(p)+=P0tuvx_no(j,k,q,l)*natorbsCI(q,p) - end do - end do - do p=1,n_act_orb - P0tuvx_no(j,k,p,l)=d(p) - end do - end do - end do - end do - ! 4th quarter - do j=1,n_act_orb - do k=1,n_act_orb - do l=1,n_act_orb - do p=1,n_act_orb - d(p)=0.D0 - end do - do p=1,n_act_orb - do q=1,n_act_orb - d(p)+=P0tuvx_no(j,k,l,q)*natorbsCI(q,p) - end do - end do - do p=1,n_act_orb - P0tuvx_no(j,k,l,p)=d(p) - end do - end do - end do - end do + call dgemm('T','N',(n_act_orb*n_act_orb*n_act_orb), n_act_orb, n_act_orb, 1.d0, & + P0tuvx, n_act_orb, natorbsCI, n_act_orb, 0.d0, & + tmp, (n_act_orb*n_act_orb*n_act_orb)) + + call dgemm('T','N',(n_act_orb*n_act_orb*n_act_orb), n_act_orb, n_act_orb, 1.d0, & + tmp, n_act_orb, natorbsCI, n_act_orb, 0.d0, & + P0tuvx_no, (n_act_orb*n_act_orb*n_act_orb)) + + call dgemm('T','N',(n_act_orb*n_act_orb*n_act_orb), n_act_orb, n_act_orb, 1.d0, & + P0tuvx_no, n_act_orb, natorbsCI, n_act_orb, 0.d0, & + tmp, (n_act_orb*n_act_orb*n_act_orb)) + + call dgemm('T','N',(n_act_orb*n_act_orb*n_act_orb), n_act_orb, n_act_orb, 1.d0, & + tmp, n_act_orb, natorbsCI, n_act_orb, 0.d0, & + P0tuvx_no, (n_act_orb*n_act_orb*n_act_orb)) + + deallocate(tmp) END_PROVIDER @@ -160,6 +103,7 @@ BEGIN_PROVIDER [real*8, one_ints_no, (mo_num,mo_num)] BEGIN_DOC ! Transformed one-e integrals END_DOC + integer :: i,j, p, q real*8 :: d(n_act_orb) one_ints_no(:,:)=mo_one_e_integrals(:,:) @@ -168,10 +112,8 @@ BEGIN_PROVIDER [real*8, one_ints_no, (mo_num,mo_num)] do j=1,mo_num do p=1,n_act_orb d(p)=0.D0 - end do - do p=1,n_act_orb do q=1,n_act_orb - d(p)+=one_ints_no(list_act(q),j)*natorbsCI(q,p) + d(p) = d(p) + one_ints_no(list_act(q),j)*natorbsCI(q,p) end do end do do p=1,n_act_orb @@ -183,8 +125,6 @@ BEGIN_PROVIDER [real*8, one_ints_no, (mo_num,mo_num)] do j=1,mo_num do p=1,n_act_orb d(p)=0.D0 - end do - do p=1,n_act_orb do q=1,n_act_orb d(p)+=one_ints_no(j,list_act(q))*natorbsCI(q,p) end do diff --git a/src/dav_general_mat/dav_general.irp.f b/src/dav_general_mat/dav_general.irp.f index a277d9ef..114476c2 100644 --- a/src/dav_general_mat/dav_general.irp.f +++ b/src/dav_general_mat/dav_general.irp.f @@ -82,7 +82,7 @@ subroutine davidson_general(u_in,H_jj,energies,dim_in,sze,N_st,N_st_diag_in,conv nproc_target = nproc double precision :: rss integer :: maxab - maxab = sze + maxab = sze m=1 disk_based = .False. @@ -204,7 +204,7 @@ subroutine davidson_general(u_in,H_jj,energies,dim_in,sze,N_st,N_st_diag_in,conv u_in(i,k) = r1*dcos(r2) enddo enddo - ! Normalize all states + ! Normalize all states do k=1,N_st_diag call normalize(u_in(:,k),sze) enddo @@ -228,20 +228,13 @@ subroutine davidson_general(u_in,H_jj,energies,dim_in,sze,N_st,N_st_diag_in,conv shift = N_st_diag*(iter-1) shift2 = N_st_diag*iter - if ((iter > 1).or.(itertot == 1)) then - ! Compute |W_k> = \sum_i |i> - ! ----------------------------------- + ! Compute |W_k> = \sum_i |i> + ! ----------------------------------- - ! Gram-Smitt to orthogonalize all new guess with the previous vectors - call ortho_qr(U,size(U,1),sze,shift2) - call ortho_qr(U,size(U,1),sze,shift2) + ! Gram-Smitt to orthogonalize all new guess with the previous vectors + call ortho_qr(U,size(U,1),sze,shift2) -! call H_S2_u_0_nstates_openmp(W(:,shift+1),U(:,shift+1),N_st_diag,sze) - call hpsi(W(:,shift+1),U(:,shift+1),N_st_diag,sze,h_mat) - else - ! Already computed in update below - continue - endif + call hpsi(W(:,shift+1),U(:,shift+1),N_st_diag,sze,h_mat) ! Compute h_kl = = ! ------------------------------------------- @@ -311,12 +304,12 @@ subroutine davidson_general(u_in,H_jj,energies,dim_in,sze,N_st,N_st_diag_in,conv do i=1,sze U(i,shift2+k) = & (lambda(k) * U(i,shift2+k) - W(i,shift2+k) ) & - /max(H_jj(i) - lambda (k),1.d-2) + /max(dabs(H_jj(i) - lambda (k)),1.d-2) * dsign(1d0,H_jj(i) - lambda (k)) enddo if (k <= N_st) then residual_norm(k) = u_dot_u(U(:,shift2+k),sze) - to_print(1,k) = lambda(k) + to_print(1,k) = lambda(k) to_print(2,k) = residual_norm(k) endif enddo @@ -324,7 +317,7 @@ subroutine davidson_general(u_in,H_jj,energies,dim_in,sze,N_st,N_st_diag_in,conv if ((itertot>1).and.(iter == 1)) then - !don't print + !don't print continue else write(*,'(1X,I3,1X,100(1X,F16.10,1X,F11.6,1X,ES11.3))') iter-1, to_print(1:2,1:N_st) @@ -333,11 +326,11 @@ subroutine davidson_general(u_in,H_jj,energies,dim_in,sze,N_st,N_st_diag_in,conv ! Check convergence if (iter > 1) then converged = dabs(maxval(residual_norm(1:N_st))) < threshold_davidson - endif - + endif + do k=1,N_st - if (residual_norm(k) > 1.e8) then + if (residual_norm(k) > 1.d8) then print *, 'Davidson failed' stop -1 endif @@ -365,13 +358,15 @@ subroutine davidson_general(u_in,H_jj,energies,dim_in,sze,N_st,N_st_diag_in,conv call dgemm('N','N', sze, N_st_diag, shift2, 1.d0, & U, size(U,1), y, size(y,1), 0.d0, u_in, size(u_in,1)) + do k=1,N_st_diag do i=1,sze U(i,k) = u_in(i,k) enddo enddo - call ortho_qr(U,size(U,1),sze,N_st_diag) - call ortho_qr(U,size(U,1),sze,N_st_diag) + + call ortho_qr(U,size(U,1),sze,N_st_diag) + do j=1,N_st_diag k=1 do while ((k