diff --git a/src/hfx_derivatives.F b/src/hfx_derivatives.F index da393480c4..733669de24 100644 --- a/src/hfx_derivatives.F +++ b/src/hfx_derivatives.F @@ -1478,9 +1478,16 @@ MODULE hfx_derivatives DO i=1,3 T2 => primitive_forces_virial(((i-1)*12+coord-1)*nsgfa(iset)*nsgfb(jset)*nsgfc(kset)*nsgfd(lset)+1:& ((i-1)*12+coord)*nsgfa(iset)*nsgfb(jset)*nsgfc(kset)*nsgfd(lset)) - CALL update_virial(nsgfa(iset), nsgfb(jset), nsgfc(kset), nsgfd(lset), & - pbd_buf, pbc_buf, pad_buf, pac_buf, fac, & - T2, tmp_virial, coord, i) + IF(with_mp2_density) THEN + CALL update_virial(nsgfa(iset), nsgfb(jset), nsgfc(kset), nsgfd(lset), & + pbd_buf, pbc_buf, pad_buf, pac_buf, fac, & + T2, tmp_virial, coord, i, & + pbd_buf_mp2, pbc_buf_mp2, pad_buf_mp2, pac_buf_mp2) + ELSE + CALL update_virial(nsgfa(iset), nsgfb(jset), nsgfc(kset), nsgfd(lset), & + pbd_buf, pbc_buf, pad_buf, pac_buf, fac, & + T2, tmp_virial, coord, i) + END IF IF( nspins == 2 ) THEN CALL update_virial(nsgfa(iset), nsgfb(jset), nsgfc(kset), nsgfd(lset), & pbd_buf_beta, pbc_buf_beta, pad_buf_beta, pac_buf_beta, fac, & @@ -2186,7 +2193,8 @@ MODULE hfx_derivatives ! ***************************************************************************** SUBROUTINE update_virial(ma_max, mb_max, mc_max, md_max, & pbd, pbc, pad, pac, fac, & - prim, tmp_virial, coord, l) + prim, tmp_virial, coord, l, & + pbd_mp2, pbc_mp2, pad_mp2, pac_mp2) INTEGER, INTENT(IN) :: ma_max, mb_max, mc_max, md_max REAL(dp), DIMENSION(*), INTENT(IN) :: pbd, pbc, pad, pac @@ -2195,26 +2203,65 @@ MODULE hfx_derivatives *md_max), INTENT(IN) :: prim REAL(dp) :: tmp_virial(3,3) INTEGER, INTENT(IN) :: coord, l + REAL(dp), DIMENSION(*), INTENT(IN), & + OPTIONAL :: pbd_mp2, pbc_mp2, pad_mp2, & + pac_mp2 INTEGER :: i, j, ma, mb, mc, md, p_index - REAL(dp) :: temp1, temp2, temp3, temp4 + LOGICAL :: with_mp2_density + REAL(dp) :: temp1, temp1_mp2, temp2, & + temp3, temp3_mp2, temp4 - p_index = 0 - temp4 = 0.0_dp - DO md = 1,md_max - DO mc = 1,mc_max - DO mb = 1,mb_max - temp1 = pbc((mc-1)*mb_max+mb) * fac - temp3 = pbd((md-1)*mb_max+mb) * fac - DO ma = 1,ma_max - p_index = p_index + 1 - temp2 = temp1 * pad((md-1)*ma_max+ma) + & - temp3 * pac((mc-1)*ma_max+ma) - temp4 = temp4 + temp2 * prim(p_index) - END DO !ma - END DO !mb - END DO !mc - END DO !md + with_mp2_density=.FALSE. + IF(PRESENT(pbd_mp2).AND.& + PRESENT(pbc_mp2).AND.& + PRESENT(pad_mp2).AND.& + PRESENT(pac_mp2)) with_mp2_density=.TRUE. + + IF(with_mp2_density) THEN + p_index = 0 + temp4 = 0.0_dp + DO md = 1,md_max + DO mc = 1,mc_max + DO mb = 1,mb_max + temp1 = pbc((mc-1)*mb_max+mb) * fac + temp3 = pbd((md-1)*mb_max+mb) * fac + temp1_mp2 = pbc_mp2((mc-1)*mb_max+mb) * fac + temp3_mp2 = pbd_mp2((md-1)*mb_max+mb) * fac + DO ma = 1,ma_max + p_index = p_index + 1 + ! HF-SCF + temp2 = temp1 * pad((md-1)*ma_max+ma) + & + temp3 * pac((mc-1)*ma_max+ma) + ! MP2+HF + temp2 = temp2 + & + pac((mc-1)*ma_max+ma) * temp3_mp2 + & + pac_mp2((mc-1)*ma_max+ma) * temp3 + & + pad((md-1)*ma_max+ma) * temp1_mp2 + & + pad_mp2((md-1)*ma_max+ma) * temp1 + temp4 = temp4 + temp2 * prim(p_index) + END DO !ma + END DO !mb + END DO !mc + END DO !md + ELSE + p_index = 0 + temp4 = 0.0_dp + DO md = 1,md_max + DO mc = 1,mc_max + DO mb = 1,mb_max + temp1 = pbc((mc-1)*mb_max+mb) * fac + temp3 = pbd((md-1)*mb_max+mb) * fac + DO ma = 1,ma_max + p_index = p_index + 1 + temp2 = temp1 * pad((md-1)*ma_max+ma) + & + temp3 * pac((mc-1)*ma_max+ma) + temp4 = temp4 + temp2 * prim(p_index) + END DO !ma + END DO !mb + END DO !mc + END DO !md + END IF j = l i = MOD(coord-1,3)+1 diff --git a/src/mp2.F b/src/mp2.F index e264bf0ef2..721b9a7472 100644 --- a/src/mp2.F +++ b/src/mp2.F @@ -231,12 +231,6 @@ MODULE mp2 CPSourceFileRef, only_ionode=.TRUE.) END IF - IF(qs_env%virial%pv_availability) THEN - CALL cp_assert(.FALSE.,cp_failure_level,cp_assertion_failed,routineP,& - "MP2 stress not implemented."//& - CPSourceFileRef, only_ionode=.TRUE.) - END IF - natom = SIZE(particle_set,1) ALLOCATE(kind_of(natom),STAT=stat) @@ -366,6 +360,7 @@ MODULE mp2 free_hfx_buffer=.TRUE. IF(calc_forces.AND.(.NOT.mp2_env%ri_mp2%free_hfx_buffer)) free_hfx_buffer=.FALSE. END IF + IF(qs_env%virial%pv_numer) free_hfx_buffer=.FALSE. IF(free_hfx_buffer) THEN CALL timeset(routineN//"_free_hfx",handle2) CALL section_vals_get(hfx_sections,n_repetition=n_rep_hf,error=error) diff --git a/src/mp2_cphf.F b/src/mp2_cphf.F index 6d7d3a810b..80d63533a0 100644 --- a/src/mp2_cphf.F +++ b/src/mp2_cphf.F @@ -54,12 +54,17 @@ MODULE mp2_cphf USE kinds, ONLY: dp USE linear_systems, ONLY: solve_system USE machine, ONLY: m_walltime + USE mathconstants, ONLY: fourpi USE message_passing, ONLY: mp_sum USE mp2_types, ONLY: mp2_type USE particle_types, ONLY: particle_type USE pw_env_types, ONLY: pw_env_get,& pw_env_type - USE pw_methods, ONLY: pw_scale,& + USE pw_methods, ONLY: pw_axpy,& + pw_copy,& + pw_derive,& + pw_integral_ab,& + pw_scale,& pw_transfer USE pw_poisson_methods, ONLY: pw_poisson_solve USE pw_poisson_types, ONLY: pw_poisson_type @@ -84,6 +89,7 @@ MODULE mp2_cphf qs_rho_type USE timings, ONLY: timeset,& timestop + USE virial_types, ONLY: virial_type #include "cp_common_uses.h" IMPLICIT NONE @@ -124,17 +130,21 @@ MODULE mp2_cphf CHARACTER(LEN=*), PARAMETER :: routineN = 'solve_z_vector_eq', & routineP = moduleN//':'//routineN - INTEGER :: bin, dimen, handle, handle2, i, i_global, i_thread, iiB, & - ikind, irep, ispin, j_global, jjB, my_bin_size, n_rep_hf, n_threads, & - ncol_local, nrow_local, stat, transf_type_in, transf_type_out, virtual + INTEGER :: alpha, beta, bin, dimen, handle, handle2, i, i_global, & + i_thread, iiB, ikind, irep, ispin, j_global, jjB, my_bin_size, & + n_rep_hf, n_threads, ncol_local, nrow_local, stat, transf_type_in, & + transf_type_out, virtual + INTEGER, DIMENSION(3) :: comp INTEGER, DIMENSION(:), POINTER :: col_indices, row_indices LOGICAL :: do_dynamic_load_balancing, & do_hfx, failure, & - hfx_treat_lsd_in_core - REAL(KIND=dp) :: out_alpha, pair_energy, & - tot_rho_r + hfx_treat_lsd_in_core, & + use_virial + REAL(KIND=dp) :: e_hartree, out_alpha, & + pair_energy, tot_rho_r REAL(KIND=dp), ALLOCATABLE, & DIMENSION(:, :) :: mat_deb + REAL(KIND=dp), DIMENSION(3, 3) :: h_stress TYPE(cp_blacs_env_type), POINTER :: blacs_env TYPE(cp_dbcsr_p_type) :: P_mu_nu TYPE(cp_dbcsr_p_type), DIMENSION(:), & @@ -148,7 +158,8 @@ MODULE mp2_cphf TYPE(hfx_container_type), POINTER :: maxval_container TYPE(hfx_type), POINTER :: actual_x_data TYPE(pw_env_type), POINTER :: pw_env - TYPE(pw_p_type) :: pot_g, rho_g, rho_r + TYPE(pw_p_type) :: dvg(3), pot_g, rho_g, rho_r, & + temp_pw_g TYPE(pw_poisson_type), POINTER :: poisson_env TYPE(pw_pool_p_type), DIMENSION(:), & POINTER :: pw_pools @@ -158,6 +169,7 @@ MODULE mp2_cphf POINTER :: force TYPE(qs_rho_type), POINTER :: rho, rho_work TYPE(section_vals_type), POINTER :: hfx_sections, input + TYPE(virial_type), POINTER :: virial !$ INTEGER :: omp_get_max_threads,omp_get_thread_num @@ -167,7 +179,7 @@ MODULE mp2_cphf ! start collecting stuff dimen=nmo virtual=dimen-homo - NULLIFY(input,pw_env,matrix_s,blacs_env,rho,energy,force) + NULLIFY(input,pw_env,matrix_s,blacs_env,rho,energy,force,virial) CALL get_qs_env(qs_env=qs_env,& pw_env=pw_env,& input=input,& @@ -176,8 +188,12 @@ MODULE mp2_cphf rho=rho,& energy=energy,& force=force,& + virial=virial,& error=error) + ! check if we have to calculate the virial + use_virial = virial%pv_availability.AND.(.NOT.virial%pv_numer) + ! mp2 matrices NULLIFY(P_mo, W_mo, L_jb) P_mo => mp2_env%ri_grad%P_mo @@ -532,6 +548,70 @@ MODULE mp2_cphf ! right contribution CALL integrate_v_rspace(v_rspace=rho_r,p=qs_env%rho%rho_ao(1),h=qs_env%matrix_ks(1),& qs_env=qs_env,calculate_forces=.TRUE.,error=error) + + IF(use_virial) THEN + ! update virial if necessery with the volume term + ! first create pw auxilliary stuff + CALL timeset(routineN//"_Virial",handle2) + NULLIFY(temp_pw_g%pw) + CALL pw_pool_create_pw(auxbas_pw_pool,temp_pw_g%pw,& + use_data=COMPLEXDATA1D,& + in_space=RECIPROCALSPACE,error=error) + DO i=1, 3 + NULLIFY(dvg(i)%pw) + CALL pw_pool_create_pw(auxbas_pw_pool,dvg(i)%pw,& + use_data=COMPLEXDATA1D,& + in_space=RECIPROCALSPACE,error=error) + END DO + + ! make a copy of the MP2 density in G space + CALL pw_copy(rho_g%pw, temp_pw_g%pw, error=error) + ! calculate MP2-like-hartree potential derivatives + DO i=1, 3 + comp=0 + comp(i)=1 + CALL pw_copy(pot_g%pw, dvg(i)%pw, error=error) + CALL pw_derive(dvg(i)%pw, comp, error=error) + END DO + + ! calculate total SCF density and potential + CALL calculate_rho_elec(matrix_p=qs_env%rho%rho_ao(1)%matrix,& + rho=rho_r,& + rho_gspace=rho_g,& + total_rho=tot_rho_r,& + qs_env=qs_env,& + soft_valid=.FALSE.,& + error=error) + ! and associated potential + CALL pw_transfer(rho_r%pw, rho_g%pw, error=error) + ! don't forget the core density + CALL pw_axpy(qs_env%rho_core%pw, rho_g%pw, error=error) + CALL pw_poisson_solve(poisson_env,rho_g%pw, pair_energy, pot_g%pw,error=error) + + ! finally update virial with the volume contribution + e_hartree=pw_integral_ab(temp_pw_g%pw, pot_g%pw, error=error) + DO alpha=1, 3 + comp=0 + comp(alpha)=1 + CALL pw_copy(pot_g%pw, rho_g%pw, error=error) + CALL pw_derive(rho_g%pw, comp, error=error) + h_stress(alpha,alpha)=-e_hartree + DO beta=alpha, 3 + h_stress(alpha,beta)=h_stress(alpha,beta) & + -2.0_dp*pw_integral_ab(rho_g%pw, dvg(beta)%pw, error=error)/fourpi + h_stress (beta,alpha)=h_stress(alpha,beta) + END DO + END DO + virial%pv_virial = virial%pv_virial + h_stress/REAL(para_env%num_pe,dp) + + ! free stuff + CALL pw_pool_give_back_pw(auxbas_pw_pool,temp_pw_g%pw,error=error) + DO i=1, 3 + CALL pw_pool_give_back_pw(auxbas_pw_pool,dvg(i)%pw,error=error) + END DO + CALL timestop(handle2) + END IF + DO ispin=1, qs_env%dft_control%nspins CALL cp_dbcsr_add(qs_env%rho%rho_ao(ispin)%matrix, qs_env%matrix_p_mp2(ispin)%matrix, 1.0_dp, 1.0_dp, error) END DO diff --git a/src/mp2_ri_grad.F b/src/mp2_ri_grad.F index ed65f0d84e..7f6df70216 100644 --- a/src/mp2_ri_grad.F +++ b/src/mp2_ri_grad.F @@ -48,6 +48,7 @@ MODULE mp2_ri_grad USE input_section_types, ONLY: section_vals_get_subs_vals,& section_vals_type USE kinds, ONLY: dp + USE mathconstants, ONLY: fourpi USE message_passing, ONLY: & mp_alltoall, mp_comm_split_direct, mp_irecv, mp_isend, & mp_request_null, mp_sendrecv, mp_sum, mp_wait, mp_waitall @@ -60,15 +61,20 @@ MODULE mp2_ri_grad USE particle_types, ONLY: particle_type USE pw_env_types, ONLY: pw_env_get,& pw_env_type - USE pw_methods, ONLY: pw_scale,& + USE pw_methods, ONLY: pw_copy,& + pw_derive,& + pw_integral_ab,& + pw_scale,& pw_transfer USE pw_poisson_methods, ONLY: pw_poisson_solve USE pw_poisson_types, ONLY: pw_poisson_type USE pw_pool_types, ONLY: pw_pool_create_pw,& pw_pool_give_back_pw,& pw_pool_type - USE pw_types, ONLY: REALDATA3D,& + USE pw_types, ONLY: COMPLEXDATA1D,& + REALDATA3D,& REALSPACE,& + RECIPROCALSPACE,& pw_p_type USE qs_collocate_density, ONLY: calculate_rho_elec,& calculate_wavefunction @@ -87,6 +93,7 @@ MODULE mp2_ri_grad USE timings, ONLY: timeset,& timestop USE util, ONLY: get_limit + USE virial_types, ONLY: virial_type #include "cp_common_uses.h" IMPLICIT NONE @@ -141,22 +148,25 @@ MODULE mp2_ri_grad CHARACTER(LEN=*), PARAMETER :: routineN = 'calc_ri_mp2_nonsep', & routineP = moduleN//':'//routineN - INTEGER :: atom_a, dimen, dir, handle, handle2, handle3, i, iatom, & - igrid_level, ikind, iorb, ipgf, iset, istat, itmp(2), L_counter, lb(3), & - LLL, location(3), my_P_end, my_P_size, my_P_start, na1, na2, natom, & - ncoa, nkind, nseta, offset, sgfa, stat, tp(3), ub(3), virtual + INTEGER :: alpha, atom_a, beta, dimen, dir, handle, handle2, handle3, i, & + iatom, igrid_level, ikind, iorb, ipgf, iset, istat, itmp(2), L_counter, & + lb(3), LLL, location(3), my_P_end, my_P_size, my_P_start, na1, na2, & + natom, ncoa, nkind, nseta, offset, sgfa, stat, tp(3), ub(3), virtual INTEGER, ALLOCATABLE, DIMENSION(:) :: atom_of_kind, kind_of + INTEGER, DIMENSION(3) :: comp INTEGER, DIMENSION(:), POINTER :: la_max, la_min, lb_max, & lb_min, npgfa, nsgfa INTEGER, DIMENSION(:, :), POINTER :: first_sgfa LOGICAL :: failure, map_it_here, & - skip_shell - REAL(KIND=dp) :: eps_filter, pair_energy, & - rab2, total_rho + skip_shell, use_virial + REAL(KIND=dp) :: e_hartree, eps_filter, & + pair_energy, rab2, total_rho REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: wf_vector REAL(KIND=dp), ALLOCATABLE, & DIMENSION(:, :) :: G_PQ_local, I_ab REAL(KIND=dp), DIMENSION(3) :: force_a, force_b, ra, rab + REAL(KIND=dp), DIMENSION(3, 3) :: h_stress, my_virial_a, & + my_virial_b REAL(KIND=dp), DIMENSION(:), POINTER :: set_radius_a REAL(KIND=dp), DIMENSION(:, :), POINTER :: I_tmp2, pab, rpgfa, sphi_a, & zeta, zetb @@ -174,7 +184,7 @@ MODULE mp2_ri_grad POINTER :: molecule_kind_set TYPE(molecule_type), DIMENSION(:), & POINTER :: molecule_set - TYPE(pw_p_type) :: psi_L + TYPE(pw_p_type) :: dvg(3), psi_L, temp_pw_g TYPE(qs_force_type), DIMENSION(:), & POINTER :: force TYPE(realspace_grid_desc_p_type), & @@ -182,6 +192,7 @@ MODULE mp2_ri_grad TYPE(realspace_grid_p_type), & DIMENSION(:), POINTER :: rs_v TYPE(section_vals_type), POINTER :: input, interp_section + TYPE(virial_type), POINTER :: virial CALL timeset(routineN,handle) failure=.FALSE. @@ -268,8 +279,27 @@ MODULE mp2_ri_grad in_space=REALSPACE,error=error_sub) ! get input section and forces - NULLIFY (force) - CALL get_qs_env(qs_env=qs_env,input=input,force=force,error=error) + NULLIFY (force, virial) + CALL get_qs_env(qs_env=qs_env,input=input,force=force,virial=virial,error=error) + + ! check if we want to calculate the virial + use_virial = virial%pv_availability.AND.(.NOT.virial%pv_numer) + + ! in case virial is required we need auxilliary pw + ! for calculate the MP2-volume contribution to the virial + ! (hartree potential derivatives) + IF(use_virial) THEN + NULLIFY(temp_pw_g%pw) + CALL pw_pool_create_pw(auxbas_pw_pool,temp_pw_g%pw,& + use_data=COMPLEXDATA1D,& + in_space=RECIPROCALSPACE,error=error_sub) + DO i=1, 3 + NULLIFY(dvg(i)%pw) + CALL pw_pool_create_pw(auxbas_pw_pool,dvg(i)%pw,& + use_data=COMPLEXDATA1D,& + in_space=RECIPROCALSPACE,error=error_sub) + END DO + END IF ! start main loop over auxiliary basis functions CALL timeset(routineN//"_loop",handle2) @@ -331,6 +361,20 @@ MODULE mp2_ri_grad CALL pw_transfer(pot_g%pw, rho_r%pw, error=error_sub) CALL pw_scale(rho_r%pw,rho_r%pw%pw_grid%dvol, error=error_sub) CALL timestop(handle3) + + IF(use_virial) THEN + ! make a copy of the density in G space + ! calculate the potential derivatives in G space + CALL timeset(routineN//"_Virial",handle3) + CALL pw_copy(rho_g%pw, temp_pw_g%pw, error=error) + DO i=1, 3 + comp=0 + comp(i)=1 + CALL pw_copy(pot_g%pw, dvg(i)%pw, error=error_sub) + CALL pw_derive(dvg(i)%pw, comp, error=error_sub) + END DO + CALL timestop(handle3) + END IF ! integrate the potential of the single gaussian and update ! 2-center forces with Gamma_PQ @@ -368,6 +412,10 @@ MODULE mp2_ri_grad force_a(:) = 0.0_dp force_b(:) = 0.0_dp + IF (use_virial) THEN + my_virial_a = 0.0_dp + my_virial_b = 0.0_dp + END IF DO iset=1, nseta ncoa = npgfa(iset)*ncoset(la_max(iset)) @@ -438,6 +486,7 @@ MODULE mp2_ri_grad eps_gvg_rspace=qs_env%dft_control%qs_control%eps_gvg_rspace,& calculate_forces=.TRUE.,& force_a=force_a, force_b=force_b, & + use_virial=use_virial,my_virial_a=my_virial_a,my_virial_b=my_virial_b,& error=error_sub) END DO @@ -457,6 +506,9 @@ MODULE mp2_ri_grad force(ikind)%rho_elec(:,atom_a) =& force(ikind)%rho_elec(:,atom_a) + force_a(:) + force_b + IF(use_virial) THEN + virial%pv_virial = virial%pv_virial + my_virial_a + my_virial_b + END IF END DO DO i=1,SIZE(rs_v) @@ -501,6 +553,43 @@ MODULE mp2_ri_grad CPPostcondition(stat==0,cp_failure_level,routineP,error,failure) CALL timestop(handle3) + IF(use_virial) THEN + ! add the volume contribution to the virial due to + ! the (P|Q) integrals, first we put the full gamme_PQ + ! pseudo wave-function into grid in order to calculate the + ! hartree potential derivatives + CALL timeset(routineN//"_Virial",handle3) + wf_vector=0.0_dp + wf_vector=-2.0_dp*G_PQ_local(:,L_counter) + CALL calculate_wavefunction(mo_coeff,1,psi_L,rho_g, & + atomic_kind_set,cell,dft_control,particle_set, & + pw_env_sub,basis_set_id=use_ri_aux_basis_set,& + external_vector=wf_vector,& + error=error_sub) + ! transfer to reciprocal space and calculate potential + rho_r%pw%cr3d = psi_L%pw%cr3d + CALL pw_transfer(rho_r%pw, rho_g%pw, error=error_sub) + CALL pw_poisson_solve(poisson_env,rho_g%pw, pair_energy, pot_g%pw,error=error_sub) + ! update virial with volume term (first calculate hartree like energy (diagonal part of the virial)) + e_hartree=0.0_dp + h_stress=0.0_dp + e_hartree=pw_integral_ab(temp_pw_g%pw, pot_g%pw, error=error_sub) + DO alpha=1, 3 + comp=0 + comp(alpha)=1 + CALL pw_copy(pot_g%pw, rho_g%pw, error=error_sub) + CALL pw_derive(rho_g%pw, comp, error=error_sub) + h_stress(alpha,alpha)=-e_hartree + DO beta=alpha, 3 + h_stress(alpha,beta)=h_stress(alpha,beta) & + -2.0_dp*pw_integral_ab(rho_g%pw, dvg(beta)%pw, error=error_sub)/fourpi + h_stress (beta,alpha)=h_stress(alpha,beta) + END DO + END DO + virial%pv_virial = virial%pv_virial + h_stress/REAL(para_env_sub%num_pe,dp) + CALL timestop(handle3) + END IF + ! put the gamma density on grid CALL timeset(routineN//"_Gpot",handle3) CALL calculate_rho_elec(matrix_p=matrix_P_munu%matrix,& @@ -518,6 +607,29 @@ MODULE mp2_ri_grad CALL pw_scale(rho_r%pw,rho_r%pw%pw_grid%dvol, error=error_sub) CALL timestop(handle3) + IF(use_virial) THEN + ! add the volume contribution to the virial coming from + ! the 3-center integrals (mu nu|P) + CALL timeset(routineN//"_Virial",handle3) + e_hartree=0.0_dp + h_stress=0.0_dp + e_hartree=pw_integral_ab(temp_pw_g%pw, pot_g%pw, error=error_sub) + DO alpha=1, 3 + comp=0 + comp(alpha)=1 + CALL pw_copy(pot_g%pw, rho_g%pw, error=error_sub) + CALL pw_derive(rho_g%pw, comp, error=error_sub) + h_stress(alpha,alpha)=-e_hartree + DO beta=alpha, 3 + h_stress(alpha,beta)=h_stress(alpha,beta) & + -2.0_dp*pw_integral_ab(rho_g%pw, dvg(beta)%pw, error=error_sub)/fourpi + h_stress (beta,alpha)=h_stress(alpha,beta) + END DO + END DO + virial%pv_virial = virial%pv_virial + h_stress/REAL(para_env_sub%num_pe,dp) + CALL timestop(handle3) + END IF + ! integrate potential with auxiliary basis function derivatives NULLIFY(rs_v) NULLIFY(rs_descs) @@ -552,6 +664,10 @@ MODULE mp2_ri_grad force_a(:) = 0.0_dp force_b(:) = 0.0_dp + IF (use_virial) THEN + my_virial_a = 0.0_dp + my_virial_b = 0.0_dp + END IF DO iset=1, nseta ncoa = npgfa(iset)*ncoset(la_max(iset)) @@ -635,6 +751,7 @@ MODULE mp2_ri_grad eps_gvg_rspace=qs_env%dft_control%qs_control%eps_gvg_rspace,& calculate_forces=.TRUE.,& force_a=force_a, force_b=force_b, & + use_virial=use_virial,my_virial_a=my_virial_a,my_virial_b=my_virial_b,& error=error_sub) END DO @@ -654,6 +771,9 @@ MODULE mp2_ri_grad force(ikind)%rho_elec(:,atom_a) =& force(ikind)%rho_elec(:,atom_a) + force_a(:) + force_b(:) + IF(use_virial) THEN + virial%pv_virial = virial%pv_virial + my_virial_a + my_virial_b + END IF END DO DO i=1,SIZE(rs_v) @@ -670,6 +790,13 @@ MODULE mp2_ri_grad CPPostcondition(stat==0,cp_failure_level,routineP,error,failure) CALL pw_pool_give_back_pw(auxbas_pw_pool,psi_L%pw,error=error_sub) + IF(use_virial) THEN + CALL pw_pool_give_back_pw(auxbas_pw_pool,temp_pw_g%pw,error=error_sub) + DO i=1, 3 + CALL pw_pool_give_back_pw(auxbas_pw_pool,dvg(i)%pw,error=error_sub) + END DO + END IF + CALL cp_dbcsr_release(matrix_P_inu%matrix,error=error_sub) DEALLOCATE(matrix_P_inu%matrix,STAT=stat) CPPostcondition(stat==0,cp_failure_level,routineP,error,failure) diff --git a/src/qs_energy.F b/src/qs_energy.F index df61c9c143..f130e77ff9 100644 --- a/src/qs_energy.F +++ b/src/qs_energy.F @@ -138,6 +138,10 @@ CONTAINS qs_env=qs_env, & calculate_forces=.FALSE., & just_energy=.TRUE.,error=error) + ! add MP2 energy if necessary + IF(ASSOCIATED(qs_env%mp2_env)) THEN + qs_env%energy%total = qs_env%energy%total + qs_env%energy%mp2 + END IF END IF END IF diff --git a/tests/QS/regtest-mp2-stress/H2O_stress_an.inp b/tests/QS/regtest-mp2-stress/H2O_stress_an.inp new file mode 100644 index 0000000000..a2ce8f4d4e --- /dev/null +++ b/tests/QS/regtest-mp2-stress/H2O_stress_an.inp @@ -0,0 +1,103 @@ +&GLOBAL + PROJECT STRESS_RI-MP2 + PRINT_LEVEL LOW + RUN_TYPE CELL_OPT +&END GLOBAL +&MOTION + &CELL_OPT + MAX_ITER 1 + &END +&END MOTION +&FORCE_EVAL + METHOD Quickstep + STRESS_TENSOR ANALYTICAL + &PRINT + &FORCES + &END + &STRESS_TENSOR + &END + &END + &DFT + BASIS_SET_FILE_NAME ../HFX_BASIS + POTENTIAL_FILE_NAME ../GTH_POTENTIALS + &MGRID + CUTOFF 150 + REL_CUTOFF 30 + &END MGRID + &QS + METHOD GPW + EPS_DEFAULT 1.0E-12 + &END QS + &SCF + SCF_GUESS ATOMIC + EPS_SCF 1.0E-6 + MAX_SCF 100 + &PRINT + &RESTART OFF + &END + &END + &END SCF + &XC + &XC_FUNCTIONAL NONE + &END XC_FUNCTIONAL + &HF + FRACTION 1.0000000 + &SCREENING + EPS_SCHWARZ 1.0E-9 + EPS_SCHWARZ_FORCES 1.0E-9 + SCREEN_ON_INITIAL_P .FALSE. + &END SCREENING + &INTERACTION_POTENTIAL + POTENTIAL_TYPE TRUNCATED + CUTOFF_RADIUS 2.4999999 + T_C_G_DATA ../t_c_g.dat + &END + &END HF + &WF_CORRELATION + METHOD RI_MP2_GPW + &RI_MP2 + BLOCK_SIZE 1 + EPS_CANONICAL 0.0001 + FREE_HFX_BUFFER .TRUE. + &END + &CPHF + EPS_CONV 1.0E-6 + MAX_ITER 50 + &END + &WFC_GPW + CUTOFF 150 + REL_CUTOFF 30 + EPS_FILTER 1.0E-12 + EPS_GRID 1.0E-8 + &END + MEMORY 1.00 + NUMBER_PROC 1 + &END + &END XC + &END DFT + &SUBSYS + &CELL + ABC [angstrom] 5.000 5.000 5.000 + &END CELL + &KIND H + BASIS_SET DZVP-GTH + RI_AUX_BASIS_SET RI_DZVP-GTH + POTENTIAL GTH-HF-q1 + &END KIND + &KIND O + BASIS_SET DZVP-GTH + RI_AUX_BASIS_SET RI_DZVP-GTH + POTENTIAL GTH-HF-q6 + &END KIND + &COORD + O 0.000000 0.000000 -0.211000 + H 0.000000 -0.844000 0.495000 + H 0.000000 0.744000 0.495000 + &END + &TOPOLOGY + &CENTER_COORDINATES + &END + &END TOPOLOGY + &END SUBSYS +&END FORCE_EVAL + diff --git a/tests/QS/regtest-mp2-stress/H2_stress_num.inp b/tests/QS/regtest-mp2-stress/H2_stress_num.inp new file mode 100644 index 0000000000..0e6d84e4a0 --- /dev/null +++ b/tests/QS/regtest-mp2-stress/H2_stress_num.inp @@ -0,0 +1,97 @@ +&GLOBAL + PROJECT STRESS_NUM_RI-MP2 + PRINT_LEVEL LOW + RUN_TYPE CELL_OPT +&END GLOBAL +&MOTION + &CELL_OPT + MAX_ITER 1 + &END +&END MOTION +&FORCE_EVAL + METHOD Quickstep + STRESS_TENSOR DIAGONAL_NUMERICAL + &PRINT + &FORCES + &END + &STRESS_TENSOR + &END + &END + &DFT + BASIS_SET_FILE_NAME ../HFX_BASIS + POTENTIAL_FILE_NAME ../GTH_POTENTIALS + &MGRID + CUTOFF 150 + REL_CUTOFF 30 + &END MGRID + &QS + METHOD GPW + EPS_DEFAULT 1.0E-12 + &END QS + &SCF + SCF_GUESS ATOMIC + EPS_SCF 1.0E-10 + MAX_SCF 100 + &PRINT + &RESTART OFF + &END + &END + &END SCF + &XC + &XC_FUNCTIONAL NONE + &END XC_FUNCTIONAL + &HF + FRACTION 1.0000000 + &SCREENING + EPS_SCHWARZ 1.0E-9 + EPS_SCHWARZ_FORCES 1.0E-9 + SCREEN_ON_INITIAL_P .FALSE. + &END SCREENING + &INTERACTION_POTENTIAL + POTENTIAL_TYPE TRUNCATED + CUTOFF_RADIUS 2.4999999 + T_C_G_DATA ../t_c_g.dat + &END + &END HF + &WF_CORRELATION + METHOD RI_MP2_GPW + &RI_MP2 + BLOCK_SIZE 1 + EPS_CANONICAL 0.0001 + FREE_HFX_BUFFER .TRUE. + &END + &CPHF + EPS_CONV 1.0E-6 + MAX_ITER 50 + &END + &WFC_GPW + CUTOFF 100 + REL_CUTOFF 20 + EPS_FILTER 1.0E-12 + EPS_GRID 1.0E-8 + &END + MEMORY 1.00 + NUMBER_PROC 1 + &END + &END XC + &END DFT + &SUBSYS + &CELL + ABC [angstrom] 5.000 5.000 5.000 + &END CELL + &KIND H + BASIS_SET DZVP-GTH + RI_AUX_BASIS_SET RI_DZVP-GTH + POTENTIAL GTH-HF-q1 + &END KIND + &COORD + H 0.000000 0.000000 0.000000 + H 0.000000 0.000000 0.740000 + &END + &TOPOLOGY + &CENTER_COORDINATES + &END + &END TOPOLOGY + &END SUBSYS +&END FORCE_EVAL + diff --git a/tests/QS/regtest-mp2-stress/TEST_FILES b/tests/QS/regtest-mp2-stress/TEST_FILES new file mode 100644 index 0000000000..ba264e0306 --- /dev/null +++ b/tests/QS/regtest-mp2-stress/TEST_FILES @@ -0,0 +1,2 @@ +H2O_stress_an.inp 11 +H2_stress_num.inp 11 diff --git a/tests/QS/regtest-mp2-stress/TEST_FILES_RESET b/tests/QS/regtest-mp2-stress/TEST_FILES_RESET new file mode 100644 index 0000000000..e69de29bb2 diff --git a/tests/TEST_DIRS b/tests/TEST_DIRS index 64a6c67b68..5355e3554c 100644 --- a/tests/TEST_DIRS +++ b/tests/TEST_DIRS @@ -2,6 +2,7 @@ # Directories have been reordered according the execution time needed for a gfortran pdbg run using 2 MPI tasks # in case a new directory is added just add it at the top of the list.. # the order will be regularly checked and modified... +QS/regtest-mp2-stress QS/regtest-mp2-grad Fist/regtest-7-1 QS/regtest-dft-vdw-corr