comp-chem/cp2k
1
0
Fork 0
mirror of https://github.com/cp2k/cp2k.git synced 2026-07-28 22:25:32 -04:00
Code Issues Projects Releases Packages Wiki Activity Actions

RI-MP2 stress tensor (Analytical and Numerical)

Browse source 
svn-origin-rev: 13497
This commit is contained in:
Mauro Del Ben 2014-01-24 14:29:47 +00:00
parent 0f8131c708
commit fc2df5b857
10 changed files with 504 additions and 48 deletions
Download patch file Download diff file Expand all files Collapse all files

89
src/hfx_derivatives.F
View file

@ -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

7
src/mp2.F
View file

@ -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)

98
src/mp2_cphf.F
View file

@ -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

151
src/mp2_ri_grad.F
View file

@ -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)

4
src/qs_energy.F
View file

@ -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

103
tests/QS/regtest-mp2-stress/H2O_stress_an.inp Normal file
View file

@ -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

97
tests/QS/regtest-mp2-stress/H2_stress_num.inp Normal file
View file

@ -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

2
tests/QS/regtest-mp2-stress/TEST_FILES Normal file
View file

@ -0,0 +1,2 @@
H2O_stress_an.inp 11
H2_stress_num.inp 11

0
tests/QS/regtest-mp2-stress/TEST_FILES_RESET Normal file
View file

1
tests/TEST_DIRS
View file

@ -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