diff --git a/src/becke_numerical_grid/example.irp.f b/src/becke_numerical_grid/example.irp.f index d91e1bc9..b9238fef 100644 --- a/src/becke_numerical_grid/example.irp.f +++ b/src/becke_numerical_grid/example.irp.f @@ -50,6 +50,7 @@ subroutine example_becke_numerical_grid print*,'The second example uses the grid points as a collection of spherical grids centered on each atom' print*,'This is mostly useful if one needs to split contributions between radial/angular/atomic of an integral' ! you browse the nuclei + integral_2 = 0.d0 do i = 1, nucl_num ! you browse the radial points attached to each nucleus do j = 1, n_points_radial_grid diff --git a/src/becke_numerical_grid/grid_becke_vector.irp.f b/src/becke_numerical_grid/grid_becke_vector.irp.f index 9da8a099..f696b023 100644 --- a/src/becke_numerical_grid/grid_becke_vector.irp.f +++ b/src/becke_numerical_grid/grid_becke_vector.irp.f @@ -71,8 +71,8 @@ END_PROVIDER enddo enddo - FREE grid_points_per_atom - FREE final_weight_at_r +! FREE grid_points_per_atom +! FREE final_weight_at_r call wall_time(wall1) print *, ' wall time for final_grid_points,', wall1 - wall0 diff --git a/src/casscf_cipsi/neworbs.irp.f b/src/casscf_cipsi/neworbs.irp.f index ca2deebb..b44ac065 100644 --- a/src/casscf_cipsi/neworbs.irp.f +++ b/src/casscf_cipsi/neworbs.irp.f @@ -116,7 +116,8 @@ END_PROVIDER end if end do if(best_vector_ovrlp_casscf.lt.0)then - best_vector_ovrlp_casscf = minloc(SXeigenval,nMonoEx+1) +! best_vector_ovrlp_casscf = minloc(SXeigenval) + best_vector_ovrlp_casscf = 1 endif c0=SXeigenvec(1,best_vector_ovrlp_casscf) if (bavard) then diff --git a/src/determinants/tr_density_matrix.irp.f b/src/determinants/tr_density_matrix.irp.f index 6faef219..856e1dc0 100644 --- a/src/determinants/tr_density_matrix.irp.f +++ b/src/determinants/tr_density_matrix.irp.f @@ -1,5 +1,5 @@ BEGIN_PROVIDER [double precision, one_e_tr_dm_mo, (mo_num, mo_num, N_states, N_states)] - + use OMP_LIB implicit none BEGIN_DOC @@ -16,6 +16,19 @@ BEGIN_PROVIDER [double precision, one_e_tr_dm_mo, (mo_num, mo_num, N_states, N_s double precision, allocatable :: tmp_a(:,:,:,:), tmp_b(:,:,:,:) integer :: krow, kcol, lrow, lcol + double precision :: mem_tot_tr_dm + integer :: nthreads + nthreads = OMP_GET_MAX_THREADS() + mem_tot_tr_dm = 8.d0*mo_num*mo_num*n_states*n_states*1d-9*(nthreads+1) + if(mem_tot_tr_dm .gt. 0.9d0 * qp_max_mem) then + print*,'Warning ! you are providing one_e_tr_dm_mo in parallel ! ' + print*,'Using that amount of threads ',nthreads + print*,'Each thread needs that amount of Gb ',8.d0*mo_num*mo_num*n_states*n_states*1d-9 + print*,'So it takes in total that amount of Gb ',mem_tot_tr_dm + print*,'The max memory for the calculation is ',qp_max_mem + print*,'It may crash because of lack of memory ...' + endif + PROVIDE psi_det_alpha_unique psi_det_beta_unique one_e_tr_dm_mo = 0d0 diff --git a/src/dft_utils_in_r/dipole_becke.irp.f b/src/dft_utils_in_r/dipole_becke.irp.f new file mode 100644 index 00000000..f82faf90 --- /dev/null +++ b/src/dft_utils_in_r/dipole_becke.irp.f @@ -0,0 +1,23 @@ + BEGIN_PROVIDER [ double precision, mo_dipole_x_becke, (mo_num,mo_num)] +&BEGIN_PROVIDER [ double precision, mo_dipole_y_becke, (mo_num,mo_num)] +&BEGIN_PROVIDER [ double precision, mo_dipole_z_becke, (mo_num,mo_num)] + implicit none + integer :: i,j,ipoint + double precision :: r(3), weight + mo_dipole_x_becke=0.d0 + mo_dipole_y_becke=0.d0 + mo_dipole_z_becke=0.d0 + do i = 1, mo_num + do j = 1, mo_num + do ipoint = 1, n_points_final_grid + r(1) = final_grid_points(1,i) + r(2) = final_grid_points(2,i) + r(3) = final_grid_points(3,i) + weight = final_weight_at_r_vector(i) + mo_dipole_x_becke(j,i) += r(1) * weight * mos_in_r_array_transp(ipoint,i) * mos_in_r_array_transp(ipoint,j) + mo_dipole_y_becke(j,i) += r(2) * weight * mos_in_r_array_transp(ipoint,i) * mos_in_r_array_transp(ipoint,j) + mo_dipole_z_becke(j,i) += r(3) * weight * mos_in_r_array_transp(ipoint,i) * mos_in_r_array_transp(ipoint,j) + enddo + enddo + enddo +END_PROVIDER diff --git a/src/iterations/print_summary.irp.f b/src/iterations/print_summary.irp.f index 8e6285e2..3d36d964 100644 --- a/src/iterations/print_summary.irp.f +++ b/src/iterations/print_summary.irp.f @@ -69,7 +69,9 @@ subroutine print_summary(e_,pt2_data,pt2_data_err,n_det_,n_configuration_,n_st,s print *, '< S^2 > = ', s2_(k) print *, 'E = ', e_(k) print *, 'Variance = ', pt2_data % variance(k), ' +/- ', pt2_data_err % variance(k) - print *, 'PT norm = ', dsqrt(pt2_data % overlap(k,k)), ' +/- ', 0.5d0*dsqrt(pt2_data % overlap(k,k)) * pt2_data_err % overlap(k,k) / (pt2_data % overlap(k,k)) + if(dabs(pt2_data % overlap(k,k)).gt.0.d0)then + print *, 'PT norm = ', dsqrt(pt2_data % overlap(k,k)), ' +/- ', 0.5d0*dsqrt(pt2_data % overlap(k,k)) * pt2_data_err % overlap(k,k) / (pt2_data % overlap(k,k)) + endif print *, 'PT2 = ', pt2_data % pt2(k), ' +/- ', pt2_data_err % pt2(k) print *, 'rPT2 = ', pt2_data % rpt2(k), ' +/- ', pt2_data_err % rpt2(k) print *, 'E+PT2 '//pt2_string//' = ', e_(k)+pt2_data % pt2(k), ' +/- ', pt2_data_err % pt2(k) diff --git a/src/mol_properties/EZFIO.cfg b/src/mol_properties/EZFIO.cfg index 3ddba227..be1c74f4 100644 --- a/src/mol_properties/EZFIO.cfg +++ b/src/mol_properties/EZFIO.cfg @@ -1,3 +1,10 @@ + +[dft_grid_mo_dipole] +type: logical +doc: If true, you compute the MO x/y/z dipole with DFT's grid +interface: ezfio,provider,ocaml +default: false + [print_all_transitions] type: logical doc: If true, print the transition between all the states diff --git a/src/mol_properties/NEED b/src/mol_properties/NEED index 8d89a452..50cf3a74 100644 --- a/src/mol_properties/NEED +++ b/src/mol_properties/NEED @@ -1,2 +1,3 @@ determinants davidson_undressed +dft_utils_in_r diff --git a/src/mol_properties/mo_dipole_properties.irp.f b/src/mol_properties/mo_dipole_properties.irp.f new file mode 100644 index 00000000..d60fb39a --- /dev/null +++ b/src/mol_properties/mo_dipole_properties.irp.f @@ -0,0 +1,14 @@ + BEGIN_PROVIDER [double precision, mo_prop_dipole_x , (mo_num,mo_num)] +&BEGIN_PROVIDER [double precision, mo_prop_dipole_y , (mo_num,mo_num)] +&BEGIN_PROVIDER [double precision, mo_prop_dipole_z , (mo_num,mo_num)] + implicit none + if(dft_grid_mo_dipole)then + mo_prop_dipole_x=mo_dipole_x_becke + mo_prop_dipole_y=mo_dipole_y_becke + mo_prop_dipole_z=mo_dipole_z_becke + else + mo_prop_dipole_x=mo_dipole_x + mo_prop_dipole_y=mo_dipole_y + mo_prop_dipole_z=mo_dipole_z + endif +END_PROVIDER diff --git a/src/mol_properties/multi_s_dipole_moment.irp.f b/src/mol_properties/multi_s_dipole_moment.irp.f index 8aae3bf4..798a620d 100644 --- a/src/mol_properties/multi_s_dipole_moment.irp.f +++ b/src/mol_properties/multi_s_dipole_moment.irp.f @@ -39,24 +39,28 @@ ! p,q: general spatial MOs ! gamma^{nm}: density matrix \bra{\Psi^n} a^{\dagger}_a a_i \ket{\Psi^m} END_DOC - + USE OMP_LIB integer :: istate, jstate ! States integer :: i, j ! general spatial MOs double precision :: nuclei_part_x, nuclei_part_y, nuclei_part_z + double precision :: mem_tot_tr_dm + integer :: nthreads + nthreads = OMP_GET_MAX_THREADS() + mem_tot_tr_dm = 8.d0*mo_num*mo_num*n_states*n_states*1d-9*(nthreads+1) multi_s_x_dipole_moment = 0.d0 multi_s_y_dipole_moment = 0.d0 multi_s_z_dipole_moment = 0.d0 - if(8.d0*mo_num*mo_num*n_states*n_states*1d-9 .lt. 200.d0) then + if(mem_tot_tr_dm .lt. 0.9d0 * qp_max_mem) then do jstate = 1, N_states do istate = 1, N_states do i = 1, mo_num do j = 1, mo_num - multi_s_x_dipole_moment(istate,jstate) -= one_e_tr_dm_mo(j,i,istate,jstate) * mo_dipole_x(j,i) - multi_s_y_dipole_moment(istate,jstate) -= one_e_tr_dm_mo(j,i,istate,jstate) * mo_dipole_y(j,i) - multi_s_z_dipole_moment(istate,jstate) -= one_e_tr_dm_mo(j,i,istate,jstate) * mo_dipole_z(j,i) + multi_s_x_dipole_moment(istate,jstate) -= one_e_tr_dm_mo(j,i,istate,jstate) * mo_prop_dipole_x(j,i) + multi_s_y_dipole_moment(istate,jstate) -= one_e_tr_dm_mo(j,i,istate,jstate) * mo_prop_dipole_y(j,i) + multi_s_z_dipole_moment(istate,jstate) -= one_e_tr_dm_mo(j,i,istate,jstate) * mo_prop_dipole_z(j,i) enddo enddo enddo @@ -66,6 +70,7 @@ ! no enouph memory ! on the fly scheme + print*,'Computing on the fly the various dipole matrix elements ' PROVIDE psi_det_alpha_unique psi_det_beta_unique @@ -85,7 +90,7 @@ !$OMP psi_bilinear_matrix_transp_rows, psi_bilinear_matrix_transp_columns, & !$OMP psi_det_alpha_unique, psi_det_beta_unique, & !$OMP psi_bilinear_matrix_values, psi_bilinear_matrix_transp_values, & - !$OMP mo_dipole_x, mo_dipole_y, mo_dipole_z, & + !$OMP mo_prop_dipole_x, mo_prop_dipole_y, mo_prop_dipole_z, & !$OMP multi_s_x_dipole_moment, multi_s_y_dipole_moment, multi_s_z_dipole_moment) !$OMP DO COLLAPSE(2) do istate = 1, N_states @@ -103,9 +108,9 @@ ck = psi_bilinear_matrix_values(k_a,istate)*psi_bilinear_matrix_values(k_a,jstate) do l = 1, elec_alpha_num j = occ(l,1) - multi_s_x_dipole_moment(istate,jstate) -= ck * mo_dipole_x(j,j) - multi_s_y_dipole_moment(istate,jstate) -= ck * mo_dipole_y(j,j) - multi_s_z_dipole_moment(istate,jstate) -= ck * mo_dipole_z(j,j) + multi_s_x_dipole_moment(istate,jstate) -= ck * mo_prop_dipole_x(j,j) + multi_s_y_dipole_moment(istate,jstate) -= ck * mo_prop_dipole_y(j,j) + multi_s_z_dipole_moment(istate,jstate) -= ck * mo_prop_dipole_z(j,j) enddo if (k_a == N_det) cycle @@ -121,13 +126,13 @@ call get_single_excitation_spin(tmp_det(1,1), tmp_det2, exc, phase, N_int) call decode_exc_spin(exc, h1, p1, h2, p2) ckl = psi_bilinear_matrix_values(k_a,istate)*psi_bilinear_matrix_values(l,jstate) * phase - multi_s_x_dipole_moment(istate,jstate) -= ckl * mo_dipole_x(h1,p1) - multi_s_y_dipole_moment(istate,jstate) -= ckl * mo_dipole_y(h1,p1) - multi_s_z_dipole_moment(istate,jstate) -= ckl * mo_dipole_z(h1,p1) + multi_s_x_dipole_moment(istate,jstate) -= ckl * mo_prop_dipole_x(h1,p1) + multi_s_y_dipole_moment(istate,jstate) -= ckl * mo_prop_dipole_y(h1,p1) + multi_s_z_dipole_moment(istate,jstate) -= ckl * mo_prop_dipole_z(h1,p1) ckl = psi_bilinear_matrix_values(k_a,jstate)*psi_bilinear_matrix_values(l,istate) * phase - multi_s_x_dipole_moment(istate,jstate) -= ckl * mo_dipole_x(p1,h1) - multi_s_y_dipole_moment(istate,jstate) -= ckl * mo_dipole_y(p1,h1) - multi_s_z_dipole_moment(istate,jstate) -= ckl * mo_dipole_z(p1,h1) + multi_s_x_dipole_moment(istate,jstate) -= ckl * mo_prop_dipole_x(p1,h1) + multi_s_y_dipole_moment(istate,jstate) -= ckl * mo_prop_dipole_y(p1,h1) + multi_s_z_dipole_moment(istate,jstate) -= ckl * mo_prop_dipole_z(p1,h1) endif l = l+1 if (l > N_det) exit @@ -148,9 +153,9 @@ ck = psi_bilinear_matrix_transp_values(k_b,istate)*psi_bilinear_matrix_transp_values(k_b,jstate) do l = 1, elec_beta_num j = occ(l,2) - multi_s_x_dipole_moment(istate,jstate) -= ck * mo_dipole_x(j,j) - multi_s_y_dipole_moment(istate,jstate) -= ck * mo_dipole_y(j,j) - multi_s_z_dipole_moment(istate,jstate) -= ck * mo_dipole_z(j,j) + multi_s_x_dipole_moment(istate,jstate) -= ck * mo_prop_dipole_x(j,j) + multi_s_y_dipole_moment(istate,jstate) -= ck * mo_prop_dipole_y(j,j) + multi_s_z_dipole_moment(istate,jstate) -= ck * mo_prop_dipole_z(j,j) enddo if (k_b == N_det) cycle @@ -166,13 +171,13 @@ call get_single_excitation_spin(tmp_det(1,2), tmp_det2, exc, phase, N_int) call decode_exc_spin(exc, h1, p1, h2, p2) ckl = psi_bilinear_matrix_transp_values(k_b,istate)*psi_bilinear_matrix_transp_values(l,jstate) * phase - multi_s_x_dipole_moment(istate,jstate) -= ckl * mo_dipole_x(h1,p1) - multi_s_y_dipole_moment(istate,jstate) -= ckl * mo_dipole_y(h1,p1) - multi_s_z_dipole_moment(istate,jstate) -= ckl * mo_dipole_z(h1,p1) + multi_s_x_dipole_moment(istate,jstate) -= ckl * mo_prop_dipole_x(h1,p1) + multi_s_y_dipole_moment(istate,jstate) -= ckl * mo_prop_dipole_y(h1,p1) + multi_s_z_dipole_moment(istate,jstate) -= ckl * mo_prop_dipole_z(h1,p1) ckl = psi_bilinear_matrix_transp_values(k_b,jstate)*psi_bilinear_matrix_transp_values(l,istate) * phase - multi_s_x_dipole_moment(istate,jstate) -= ckl * mo_dipole_x(p1,h1) - multi_s_y_dipole_moment(istate,jstate) -= ckl * mo_dipole_y(p1,h1) - multi_s_z_dipole_moment(istate,jstate) -= ckl * mo_dipole_z(p1,h1) + multi_s_x_dipole_moment(istate,jstate) -= ckl * mo_prop_dipole_x(p1,h1) + multi_s_y_dipole_moment(istate,jstate) -= ckl * mo_prop_dipole_y(p1,h1) + multi_s_z_dipole_moment(istate,jstate) -= ckl * mo_prop_dipole_z(p1,h1) endif l = l+1 if (l > N_det) exit