Merge branch 'dev-stable' of github.com:QuantumPackage/qp2 into dev-s…

…table

plugins/local/spher_harm/spher_harm.irp.f

@@ -1,7 +1,7 @@
 program spher_harm
  implicit none
-! call test_spher_harm
+ call test_spher_harm
 ! call test_cart
- call test_brutal_spheric
+! call test_brutal_spheric
 end
 

plugins/local/spher_harm/spher_harm_func.irp.f

@@ -7,6 +7,7 @@ subroutine spher_harm_func_r3(r,l,m,re_ylm, im_ylm)
  double precision :: theta, phi,r_abs
  call cartesian_to_spherical(r,theta,phi,r_abs)
  call spher_harm_func(l,m,theta,phi,re_ylm, im_ylm)
+! call spher_harm_func_expl(l,m,theta,phi,re_ylm, im_ylm)
 end
 
 
@@ -131,6 +132,10 @@ subroutine spher_harm_func_expl(l,m,theta,phi,re_ylm, im_ylm)
  tmp = - inv_sq_pi * dsqrt(3.d0/8.d0) * dsin(theta) 
  re_ylm = tmp * dcos(phi)
  im_ylm = tmp * dsin(phi)
+ else if (l==1.and.m==-1)then
+ tmp = - inv_sq_pi * dsqrt(3.d0/8.d0) * dsin(theta) 
+ re_ylm = tmp * dcos(phi)
+ im_ylm = -tmp * dsin(phi)
  else if(l==1.and.m==0)then
  tmp = inv_sq_pi * dsqrt(3.d0/4.d0) * dcos(theta) 
  re_ylm = tmp 
@@ -139,10 +144,18 @@ subroutine spher_harm_func_expl(l,m,theta,phi,re_ylm, im_ylm)
  tmp = 0.25d0 * inv_sq_pi * dsqrt(0.5d0*15.d0) * dsin(theta)*dsin(theta)
  re_ylm = tmp * dcos(2.d0*phi)
  im_ylm = tmp * dsin(2.d0*phi)
+ else if(l==2.and.m==-2)then
+ tmp = 0.25d0 * inv_sq_pi * dsqrt(0.5d0*15.d0) * dsin(theta)*dsin(theta)
+ re_ylm = tmp * dcos(2.d0*phi)
+ im_ylm =-tmp * dsin(2.d0*phi)
  else if(l==2.and.m==1)then
  tmp = - inv_sq_pi * dsqrt(15.d0/8.d0) * dsin(theta) * dcos(theta)
  re_ylm = tmp * dcos(phi)
  im_ylm = tmp * dsin(phi)
+ else if(l==2.and.m==-1)then
+ tmp = - inv_sq_pi * dsqrt(15.d0/8.d0) * dsin(theta) * dcos(theta)
+ re_ylm = tmp * dcos(phi)
+ im_ylm =-tmp * dsin(phi)
  else if(l==2.and.m==0)then
  tmp = dsqrt(5.d0/4.d0) * inv_sq_pi* (1.5d0*dcos(theta)*dcos(theta)-0.5d0)
  re_ylm = tmp

src/casscf_cipsi/bielec.irp.f

@@ -1,33 +1,40 @@
-BEGIN_PROVIDER [real*8, bielec_PQxx, (mo_num, mo_num,n_core_inact_act_orb,n_core_inact_act_orb)]
+BEGIN_PROVIDER [real*8, bielec_PQxx_array, (mo_num, mo_num,n_core_inact_act_orb,n_core_inact_act_orb)]
  BEGIN_DOC
- ! bielec_PQxx : integral (pq|xx) with p,q arbitrary, x core or active
+ ! WARNING !!! Old version !!! NOT USED ANYMORE IN THE PROGRAM !!! TOO BIG TO BE STORED ON LARGE SYSTEMS !!! 
+ ! 
+ ! Replaced by the Cholesky-based function bielec_PQxx
+ !
+ ! bielec_PQxx_array : integral (pq|xx) with p,q arbitrary, x core or active
  ! indices are unshifted orbital numbers
  END_DOC
  implicit none
  integer :: i,j,ii,jj,p,q,i3,j3,t3,v3
  real*8 :: mo_two_e_integral
+ print*,''
+ print*,'Providing bielec_PQxx_array, WARNING IT CAN BE A VERY BIG ARRAY WHEN MO_NUM IS LARGE !!!'
+ print*,''
 
- bielec_PQxx(:,:,:,:) = 0.d0
+ bielec_PQxx_array(:,:,:,:) = 0.d0
  PROVIDE mo_two_e_integrals_in_map
 
  !$OMP PARALLEL DEFAULT(NONE) &
  !$OMP PRIVATE(i,ii,j,jj,i3,j3) &
- !$OMP SHARED(n_core_inact_orb,list_core_inact,mo_num,bielec_PQxx, &
+ !$OMP SHARED(n_core_inact_orb,list_core_inact,mo_num,bielec_PQxx_array, &
  !$OMP n_act_orb,mo_integrals_map,list_act)
 
  !$OMP DO
  do i=1,n_core_inact_orb
  ii=list_core_inact(i)
  do j=i,n_core_inact_orb
  jj=list_core_inact(j)
- call get_mo_two_e_integrals_i1j1(ii,jj,mo_num,bielec_PQxx(1,1,i,j),mo_integrals_map)
- bielec_PQxx(:,:,j,i)=bielec_PQxx(:,:,i,j)
+ call get_mo_two_e_integrals_i1j1(ii,jj,mo_num,bielec_PQxx_array(1,1,i,j),mo_integrals_map)
+ bielec_PQxx_array(:,:,j,i)=bielec_PQxx_array(:,:,i,j)
  end do
  do j=1,n_act_orb
  jj=list_act(j)
  j3=j+n_core_inact_orb
- call get_mo_two_e_integrals_i1j1(ii,jj,mo_num,bielec_PQxx(1,1,i,j3),mo_integrals_map)
- bielec_PQxx(:,:,j3,i)=bielec_PQxx(:,:,i,j3)
+ call get_mo_two_e_integrals_i1j1(ii,jj,mo_num,bielec_PQxx_array(1,1,i,j3),mo_integrals_map)
+ bielec_PQxx_array(:,:,j3,i)=bielec_PQxx_array(:,:,i,j3)
  end do
  end do
  !$OMP END DO
@@ -40,8 +47,8 @@
  do j=i,n_act_orb
  jj=list_act(j)
  j3=j+n_core_inact_orb
- call get_mo_two_e_integrals_i1j1(ii,jj,mo_num,bielec_PQxx(1,1,i3,j3),mo_integrals_map)
- bielec_PQxx(:,:,j3,i3)=bielec_PQxx(:,:,i3,j3)
+ call get_mo_two_e_integrals_i1j1(ii,jj,mo_num,bielec_PQxx_array(1,1,i3,j3),mo_integrals_map)
+ bielec_PQxx_array(:,:,j3,i3)=bielec_PQxx_array(:,:,i3,j3)
  end do
  end do
  !$OMP END DO
@@ -52,22 +59,29 @@
 
 
 
-BEGIN_PROVIDER [real*8, bielec_PxxQ, (mo_num,n_core_inact_act_orb,n_core_inact_act_orb, mo_num)]
+BEGIN_PROVIDER [real*8, bielec_PxxQ_array, (mo_num,n_core_inact_act_orb,n_core_inact_act_orb, mo_num)]
  BEGIN_DOC
- ! bielec_PxxQ : integral (px|xq) with p,q arbitrary, x core or active
+ ! WARNING !!! Old version !!! NOT USED ANYMORE IN THE PROGRAM !!! TOO BIG TO BE STORED ON LARGE SYSTEMS !!! 
+ ! 
+ ! Replaced by the Cholesky-based function bielec_PxxQ
+ !
+ ! bielec_PxxQ_array : integral (px|xq) with p,q arbitrary, x core or active
  ! indices are unshifted orbital numbers
  END_DOC
  implicit none
  integer :: i,j,ii,jj,p,q,i3,j3,t3,v3
  double precision, allocatable :: integrals_array(:,:)
  real*8 :: mo_two_e_integral
 
+ print*,''
+ print*,'Providing bielec_PxxQ_array, WARNING IT CAN BE A VERY BIG ARRAY WHEN MO_NUM IS LARGE !!!'
+ print*,''
  PROVIDE mo_two_e_integrals_in_map
- bielec_PxxQ = 0.d0
+ bielec_PxxQ_array = 0.d0
 
  !$OMP PARALLEL DEFAULT(NONE) &
  !$OMP PRIVATE(i,ii,j,jj,i3,j3,integrals_array) &
- !$OMP SHARED(n_core_inact_orb,list_core_inact,mo_num,bielec_PxxQ, &
+ !$OMP SHARED(n_core_inact_orb,list_core_inact,mo_num,bielec_PxxQ_array, &
  !$OMP n_act_orb,mo_integrals_map,list_act)
 
  allocate(integrals_array(mo_num,mo_num))
@@ -80,8 +94,8 @@
  call get_mo_two_e_integrals_ij(ii,jj,mo_num,integrals_array,mo_integrals_map)
  do q=1,mo_num
  do p=1,mo_num
- bielec_PxxQ(p,i,j,q)=integrals_array(p,q)
- bielec_PxxQ(p,j,i,q)=integrals_array(q,p)
+ bielec_PxxQ_array(p,i,j,q)=integrals_array(p,q)
+ bielec_PxxQ_array(p,j,i,q)=integrals_array(q,p)
  end do
  end do
  end do
@@ -91,8 +105,8 @@
  call get_mo_two_e_integrals_ij(ii,jj,mo_num,integrals_array,mo_integrals_map)
  do q=1,mo_num
  do p=1,mo_num
- bielec_PxxQ(p,i,j3,q)=integrals_array(p,q)
- bielec_PxxQ(p,j3,i,q)=integrals_array(q,p)
+ bielec_PxxQ_array(p,i,j3,q)=integrals_array(p,q)
+ bielec_PxxQ_array(p,j3,i,q)=integrals_array(q,p)
  end do
  end do
  end do
@@ -111,8 +125,8 @@
  call get_mo_two_e_integrals_ij(ii,jj,mo_num,integrals_array,mo_integrals_map)
  do q=1,mo_num
  do p=1,mo_num
- bielec_PxxQ(p,i3,j3,q)=integrals_array(p,q)
- bielec_PxxQ(p,j3,i3,q)=integrals_array(q,p)
+ bielec_PxxQ_array(p,i3,j3,q)=integrals_array(p,q)
+ bielec_PxxQ_array(p,j3,i3,q)=integrals_array(q,p)
  end do
  end do
  end do
@@ -129,10 +143,15 @@
  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
+ ! 
+ ! This array can be stored anyway. Ex: 50 active orbitals, 1500 MOs ==> 8x50^3x1500 = 1.5 Gb
  END_DOC
  implicit none
  integer :: i,j,k,p,t,u,v
  double precision, external :: mo_two_e_integral
+ double precision :: wall0, wall1 
+ call wall_time(wall0)
+ print*,'Providing bielecCI'
  PROVIDE mo_two_e_integrals_in_map
 
  !$OMP PARALLEL DO DEFAULT(NONE) &
@@ -151,5 +170,7 @@
  end do
  end do
  !$OMP END PARALLEL DO
+ call wall_time(wall1)
+ print*,'Time to provide bielecCI = ',wall1 - wall0
 
 END_PROVIDER