comp-chem/cp2k
1
0
Fork 0
mirror of https://github.com/cp2k/cp2k.git synced 2026-07-21 14:35:15 -04:00
Code Issues Projects Releases Packages Wiki Activity Actions

QMM: Add OpenMP to GEEP subroutines

Browse source 
OpenMP has been added into:
- qmmm_forces_gaussian_low
- qmmm_forces_with_gaussian_LG
- qmmm_forces_with_gaussian_LR
- qmmm_elec_with_gaussian_LG
- qmmm_elec_with_gaussian_LR
This commit is contained in:
holly-t 2020-06-05 14:01:11 +01:00 • committed by GitHub
parent 5ba09c0f20
commit 6e0731f2f8
No known key found for this signature in database
GPG key ID: 4AEE18F83AFDEB23
2 changed files with 174 additions and 193 deletions
Download patch file Download diff file Expand all files Collapse all files

27
src/qmmm_gpw_energy.F
View file

@ -582,9 +582,20 @@ CONTAINS
dr2 = pw%pw_grid%dr(2)
dr3 = pw%pw_grid%dr(3)
grid => pw%cr3d(:, :, :)
!$OMP PARALLEL DO DEFAULT(NONE) &
!$OMP SHARED(bo, gbo, grid, grid2, pw, npts, per_pot, mm_atom_index) &
!$OMP SHARED(dr1, dr2, dr3, dr1c, dr2c, dr3c, par_scheme, mm_charges, mm_particles) &
!$OMP SHARED(mm_cell, dOmmOqm, shells, para_env, IRadTyp, qmmm_spherical_cutoff) &
!$OMP PRIVATE(Imm, LIndMM, IndMM, qt, sph_chrg_factor, ra, myind) &
!$OMP PRIVATE(rt1, rt2, rt3, k, vec, ivec, xd1, xd2, xd3, ik1, ik2, ik3, ik4) &
!$OMP PRIVATE(ij1, ij2, ij3, ij4, ii1, ii2, ii3, ii4, my_k, my_j, xs1, xs2, xs3) &
!$OMP PRIVATE(p1, p2, p3, q1, q2, q3, r1, r2, r3, v1, v2, v3, v4, e1, e2, e3) &
!$OMP PRIVATE(f1, f2, f3, g1, g2, g3, h1, h2, h3, s1, s2, s3, s4, a1, a2, a3) &
!$OMP PRIVATE(b1, b2, b3, c1, c2, c3, d1, d2, d3, t1, t2, t3, t4, u1, u2, u3, val) &
!$OMP PRIVATE(rv1, rv2, rv3, abc_X, abc_X_Y)
Atoms: DO Imm = 1, SIZE(per_pot%mm_atom_index)
IF (par_scheme == do_par_atom) THEN
myind = myind + 1
myind = Imm + (IRadTyp - 1)*SIZE(per_pot%mm_atom_index)
IF (MOD(myind, para_env%num_pe) /= para_env%mepos) CYCLE
END IF
LIndMM = per_pot%mm_atom_index(Imm)
@ -712,14 +723,16 @@ CONTAINS
abc_X_Y(4) = abc_X(1, 4)*t1 + abc_X(2, 4)*t2 + abc_X(3, 4)*t3 + abc_X(4, 4)*t4
val = abc_X_Y(1)*s1 + abc_X_Y(2)*s2 + abc_X_Y(3)*s3 + abc_X_Y(4)*s4
!$OMP ATOMIC
grid2(i, j, k) = grid2(i, j, k) - val*qt
!$OMP END ATOMIC
xs1 = xs1 + dr1c
END DO
xs2 = xs2 + dr2c
END DO
END DO LoopOnGrid
END DO Atoms
!$OMP END PARALLEL DO
END DO Radius
CALL timestop(handle)
END SUBROUTINE qmmm_elec_with_gaussian_LG
@ -790,9 +803,14 @@ CONTAINS
pot => potentials(IRadTyp)%pot
dx = Pot%dx
pot0_2 => Pot%pot0_2
!$OMP PARALLEL DO DEFAULT(NONE) &
!$OMP SHARED(pot, par_scheme, para_env, mm_atom_index, mm_particles, dOmmOqm, mm_cell, qmmm_spherical_cutoff) &
!$OMP SHARED(bo, gbo, dr1, dr2, dr3, grid2, shells, pot0_2, dx, mm_charges, IRadTyp) &
!$OMP PRIVATE(myind, Imm, LIndMM, IndMM, ra, qt, sph_chrg_factor, rt1, rt2, rt3, my_k, my_j) &
!$OMP PRIVATE(rv1, rv2, rv3, rx2, rx3, r, r2, rx, Term, xs1, xs2, xs3, i, j, k, ix)
Atoms: DO Imm = 1, SIZE(pot%mm_atom_index)
IF (par_scheme == do_par_atom) THEN
myind = myind + 1
myind = Imm + (IRadTyp - 1)*SIZE(pot%mm_atom_index)
IF (MOD(myind, para_env%num_pe) /= para_env%mepos) CYCLE
END IF
LIndMM = pot%mm_atom_index(Imm)
@ -831,13 +849,16 @@ CONTAINS
+ pot0_2(2, ix)*(rx - 2.0_dp*rx2 + rx3) &
+ pot0_2(1, ix + 1)*(3.0_dp*rx2 - 2.0_dp*rx3) &
+ pot0_2(2, ix + 1)*(-rx2 + rx3)
!$OMP ATOMIC
grid2(i, j, k) = grid2(i, j, k) - Term*qt
!$OMP END ATOMIC
xs1 = xs1 + dr1
END DO
xs2 = xs2 + dr2
END DO
END DO LoopOnGrid
END DO Atoms
!$OMP END PARALLEL DO
END DO Radius
CALL timestop(handle)
END SUBROUTINE qmmm_elec_with_gaussian_LR

340
src/qmmm_gpw_forces.F
View file

@ -10,71 +10,71 @@
!> \author Teodoro Laino
! **************************************************************************************************
MODULE qmmm_gpw_forces
USE cell_types, ONLY: cell_type, &
pbc
USE cp_control_types, ONLY: dft_control_type
USE cp_log_handling, ONLY: cp_get_default_logger, &
cp_logger_type
USE cp_output_handling, ONLY: cp_print_key_finished_output, &
cp_print_key_unit_nr
USE cp_para_types, ONLY: cp_para_env_type
USE cp_spline_utils, ONLY: pw_restrict_s3, &
spline3_nopbc_interp, &
spline3_pbc_interp
USE cube_utils, ONLY: cube_info_type
USE input_constants, ONLY: do_par_atom, &
do_qmmm_coulomb, &
do_qmmm_gauss, &
do_qmmm_none, &
do_qmmm_pcharge, &
do_qmmm_swave
USE input_section_types, ONLY: section_vals_get_subs_vals, &
section_vals_type, &
section_vals_val_get
USE kinds, ONLY: dp
USE message_passing, ONLY: mp_irecv, &
mp_isend, &
mp_sum, &
mp_wait
USE mm_collocate_potential, ONLY: collocate_gf_rspace_NoPBC, &
integrate_gf_rspace_NoPBC
USE particle_types, ONLY: particle_type
USE pw_env_types, ONLY: pw_env_get, &
pw_env_type
USE pw_methods, ONLY: pw_axpy, &
pw_integral_ab, &
pw_transfer, &
pw_zero
USE pw_pool_types, ONLY: pw_pool_create_pw, &
pw_pool_give_back_pw, &
pw_pool_p_type, &
pw_pool_type, &
pw_pools_create_pws, &
pw_pools_give_back_pws
USE pw_types, ONLY: REALDATA3D, &
REALSPACE, &
pw_p_type, &
pw_type
USE qmmm_gaussian_types, ONLY: qmmm_gaussian_p_type, &
qmmm_gaussian_type
USE qmmm_gpw_energy, ONLY: qmmm_elec_with_gaussian, &
qmmm_elec_with_gaussian_LG, &
qmmm_elec_with_gaussian_LR
USE qmmm_se_forces, ONLY: deriv_se_qmmm_matrix
USE qmmm_types_low, ONLY: qmmm_env_qm_type, &
qmmm_per_pot_p_type, &
qmmm_per_pot_type, &
qmmm_pot_p_type, &
qmmm_pot_type
USE qmmm_util, ONLY: spherical_cutoff_factor
USE qmmm_tb_methods, ONLY: deriv_tb_qmmm_matrix, &
deriv_tb_qmmm_matrix_pc
USE qs_energy_types, ONLY: qs_energy_type
USE qs_environment_types, ONLY: get_qs_env, &
qs_environment_type
USE qs_ks_qmmm_types, ONLY: qs_ks_qmmm_env_type
USE qs_rho_types, ONLY: qs_rho_get, &
qs_rho_type
USE cell_types, ONLY: cell_type,&
pbc
USE cp_control_types, ONLY: dft_control_type
USE cp_log_handling, ONLY: cp_get_default_logger,&
cp_logger_type
USE cp_output_handling, ONLY: cp_print_key_finished_output,&
cp_print_key_unit_nr
USE cp_para_types, ONLY: cp_para_env_type
USE cp_spline_utils, ONLY: pw_restrict_s3,&
spline3_nopbc_interp,&
spline3_pbc_interp
USE cube_utils, ONLY: cube_info_type
USE input_constants, ONLY: do_par_atom,&
do_qmmm_coulomb,&
do_qmmm_gauss,&
do_qmmm_none,&
do_qmmm_pcharge,&
do_qmmm_swave
USE input_section_types, ONLY: section_vals_get_subs_vals,&
section_vals_type,&
section_vals_val_get
USE kinds, ONLY: dp
USE message_passing, ONLY: mp_irecv,&
mp_isend,&
mp_sum,&
mp_wait
USE mm_collocate_potential, ONLY: collocate_gf_rspace_NoPBC,&
integrate_gf_rspace_NoPBC
USE particle_types, ONLY: particle_type
USE pw_env_types, ONLY: pw_env_get,&
pw_env_type
USE pw_methods, ONLY: pw_axpy,&
pw_integral_ab,&
pw_transfer,&
pw_zero
USE pw_pool_types, ONLY: pw_pool_create_pw,&
pw_pool_give_back_pw,&
pw_pool_p_type,&
pw_pool_type,&
pw_pools_create_pws,&
pw_pools_give_back_pws
USE pw_types, ONLY: REALDATA3D,&
REALSPACE,&
pw_p_type,&
pw_type
USE qmmm_gaussian_types, ONLY: qmmm_gaussian_p_type,&
qmmm_gaussian_type
USE qmmm_gpw_energy, ONLY: qmmm_elec_with_gaussian,&
qmmm_elec_with_gaussian_LG,&
qmmm_elec_with_gaussian_LR
USE qmmm_se_forces, ONLY: deriv_se_qmmm_matrix
USE qmmm_tb_methods, ONLY: deriv_tb_qmmm_matrix,&
deriv_tb_qmmm_matrix_pc
USE qmmm_types_low, ONLY: qmmm_env_qm_type,&
qmmm_per_pot_p_type,&
qmmm_per_pot_type,&
qmmm_pot_p_type,&
qmmm_pot_type
USE qmmm_util, ONLY: spherical_cutoff_factor
USE qs_energy_types, ONLY: qs_energy_type
USE qs_environment_types, ONLY: get_qs_env,&
qs_environment_type
USE qs_ks_qmmm_types, ONLY: qs_ks_qmmm_env_type
USE qs_rho_types, ONLY: qs_rho_get,&
qs_rho_type
#include "./base/base_uses.f90"
IMPLICIT NONE
@ -333,7 +333,7 @@ CONTAINS
IF (qmmm_env%image_charge) THEN
DO iatom = 1, qmmm_env%num_image_mm_atoms
image_IndMM = qmmm_env%image_charge_pot%image_mm_list(iatom)
IF (image_IndMM .eq. IndMM) THEN
IF (image_IndMM .EQ. IndMM) THEN
Forces(:, Imm) = Forces(:, Imm) &
+ qmmm_env%image_charge_pot%image_forcesMM(:, iatom)
ENDIF
@ -415,7 +415,7 @@ CONTAINS
TYPE(cell_type), POINTER :: mm_cell
CHARACTER(len=*), PARAMETER :: routineN = 'qmmm_forces_with_gaussian', &
routineP = moduleN//':'//routineN
routineP = moduleN//':'//routineN
INTEGER :: group, handle, i, igrid, j, k, &
kind_interp, me, ngrids, request, stat
@ -617,7 +617,7 @@ CONTAINS
LOGICAL, INTENT(in) :: shells
CHARACTER(len=*), PARAMETER :: routineN = 'qmmm_force_with_gaussian_low', &
routineNb = 'qmmm_forces_gaussian_low', routineP = moduleN//':'//routineN
routineNb = 'qmmm_forces_gaussian_low', routineP = moduleN//':'//routineN
INTEGER :: handle, handle2, IGauss, ilevel, Imm, &
IndMM, IRadTyp, LIndMM, myind, &
@ -652,9 +652,16 @@ CONTAINS
ALLOCATE (xdat(2, bo(1, 1):bo(2, 1)))
ALLOCATE (ydat(2, bo(1, 2):bo(2, 2)))
ALLOCATE (zdat(2, bo(1, 3):bo(2, 3)))
!$OMP PARALLEL DO DEFAULT(NONE) &
!$OMP SHARED(pot, par_scheme, dvol, alpha, para_env, mm_atom_index, shells) &
!$OMP SHARED(mm_particles, dOmmOqm, mm_cell, height, mm_charges, qmmm_spherical_cutoff) &
!$OMP SHARED(grids, cube_info, bo, n_rep_real, eps_mm_rspace, Forces, ilevel) &
!$OMP SHARED(IGauss, pgf, IRadTyp, iw, aug_pools, auxbas_grid) &
!$OMP PRIVATE(xdat, ydat, zdat) &
!$OMP PRIVATE(Imm, LIndMM, IndMM, ra, W, force, sph_chrg_factor, myind)
Atoms: DO Imm = 1, SIZE(pot%mm_atom_index)
IF (par_scheme == do_par_atom) THEN
myind = myind + 1
myind = Imm + (IGauss - 1)*SIZE(pot%mm_atom_index) + (IRadTyp - 1)*pgf%Number_of_Gaussians
IF (MOD(myind, para_env%num_pe) /= para_env%mepos) CYCLE Atoms
END IF
LIndMM = pot%mm_atom_index(Imm)
@ -707,6 +714,7 @@ CONTAINS
iw=iw)
END IF
END DO Atoms
!$OMP END PARALLEL DO
DEALLOCATE (xdat)
DEALLOCATE (ydat)
DEALLOCATE (zdat)
@ -801,17 +809,17 @@ CONTAINS
LOGICAL :: shells
CHARACTER(len=*), PARAMETER :: routineN = 'qmmm_forces_with_gaussian_LG', &
routineP = moduleN//':'//routineN
routineP = moduleN//':'//routineN
INTEGER :: handle, i, ii1, ii2, ii3, ii4, ij1, ij2, ij3, ij4, ik1, ik2, ik3, ik4, Imm, &
IndMM, IRadTyp, ivec(3), j, k, LIndMM, my_i, my_j, my_k, myind, npts(3)
IndMM, IRadTyp, ivec(3), j, k, LIndMM, my_i, my_j, my_k, myind, npts(3)
INTEGER, DIMENSION(2, 3) :: bo, gbo
REAL(KIND=dp) :: a1, a2, a3, abc_X(4, 4), abc_X_Y(4), b1, b2, b3, c1, c2, c3, d1, d2, d3, &
dr1, dr1c, dr1i, dr2, dr2c, dr2i, dr3, dr3c, dr3i, dvol, e1, e2, e3, f1, f2, f3, fac, &
ft1, ft2, ft3, g1, g2, g3, h1, h2, h3, p1, p2, p3, q1, q2, q3, qt, r1, r2, r3, rt1, rt2, &
rt3, rv1, rv2, rv3, s1, s1d, s1o, s2, s2d, s2o, s3, s3d, s3o, s4, s4d, s4o, &
sph_chrg_factor, t1, t1d, t1o, t2, t2d, t2o, t3, t3d, t3o, t4, t4d, t4o, u1, u2, u3, v1, &
v1d, v1o, v2, v2d, v2o, v3, v3d, v3o, v4, v4d, v4o, xd1, xd2, xd3, xs1, xs2, xs3
dr1, dr1c, dr1i, dr2, dr2c, dr2i, dr3, dr3c, dr3i, dvol, e1, e2, e3, f1, f2, f3, fac, &
ft1, ft2, ft3, g1, g2, g3, h1, h2, h3, p1, p2, p3, q1, q2, q3, qt, r1, r2, r3, rt1, rt2, &
rt3, rv1, rv2, rv3, s1, s1d, s1o, s2, s2d, s2o, s3, s3d, s3o, s4, s4d, s4o, &
sph_chrg_factor, t1, t1d, t1o, t2, t2d, t2o, t3, t3d, t3o, t4, t4d, t4o, u1, u2, u3, v1, &
v1d, v1o, v2, v2d, v2o, v3, v3d, v3o, v4, v4d, v4o, xd1, xd2, xd3, xs1, xs2, xs3
REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: LForces
REAL(KIND=dp), DIMENSION(3) :: ra, val, vec
REAL(KIND=dp), DIMENSION(:, :, :), POINTER :: grid, grid2
@ -841,9 +849,24 @@ CONTAINS
dr1i = 1.0_dp/dr1
dr2i = 1.0_dp/dr2
dr3i = 1.0_dp/dr3
!$OMP PARALLEL DO DEFAULT(NONE) &
!$OMP SHARED(bo, grid, grid2, pw, npts, gbo, per_pot, mm_atom_index) &
!$OMP SHARED(dr1, dr2, dr3, dr1i, dr2i, dr3i, dr1c, dr2c, dr3c, par_scheme, mm_charges) &
!$OMP SHARED(mm_cell, dOmmOqm, dvol, shells, para_env, IRadTyp) &
!$OMP SHARED(qmmm_spherical_cutoff, mm_particles, Forces, LForces) &
!$OMP PRIVATE(qt, Imm, LIndMM, IndMM, sph_chrg_factor, ra, myind) &
!$OMP PRIVATE(rt1, rt2, rt3, ft1, ft2, ft3, my_k, my_j, my_i, xs3, xs2, xs1) &
!$OMP PRIVATE(rv3, rv2, rv1, vec, ivec, ik1, ik2, ik3, ik4, xd3, xd2, xd1) &
!$OMP PRIVATE(p1, p2, p3, q1, q2, q3, r1, r2, r3, u1, u2, u3, v1o, v2o, v3o, v4o) &
!$OMP PRIVATE(v1d, v2d, v3d, v4d, ij1, ij2, ij3, ij4, e1, e2, e3, f1, f2, f3) &
!$OMP PRIVATE(g1, g2, g3, h1, h2, h3, s1o, s2o, s3o, s4o, s1d, s2d, s3d, s4d) &
!$OMP PRIVATE(ii1, ii2, ii3, ii4, a1, a2, a3, b1, b2, b3, c1, c2, c3, d1, d2, d3) &
!$OMP PRIVATE(t1o, t2o, t3o, t4o, t1d, t2d, t3d, t4d, t1, t2, t3, t4, s1, s2, s3, s4) &
!$OMP PRIVATE(v1, v2, v3, v4, abc_x, abc_x_y, val, fac)
Atoms: DO Imm = 1, SIZE(per_pot%mm_atom_index)
IF (par_scheme == do_par_atom) THEN
myind = myind + 1
myind = Imm + (IRadTyp - 1)*SIZE(per_pot%mm_atom_index)
IF (MOD(myind, para_env%num_pe) /= para_env%mepos) CYCLE
END IF
LIndMM = per_pot%mm_atom_index(Imm)
@ -960,60 +983,6 @@ CONTAINS
t3d = -2.0_dp + 4.0_dp*c1 - 1.5_dp*c2
t4d = 0.5_dp*d2
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!
! # v then t then s
!
! for ii in 1 2 3; do
! if [[ $ii == 1 ]]; then ld=t; fi
! if [[ $ii == 2 ]]; then ld=s; fi
! if [[ $ii == 3 ]]; then ld=v; fi
! #
! for l in t s v; do
! for i in 1 2 3 4; do
! if [[ $ld == $l ]]; then
! echo "$l$i = $l${i}d*dr${ii}i"
! else
! echo "$l$i = $l${i}o"
! fi
! done
! done
! #
! for i in 1 2 3 4; do
! for j in 1 2 3 4; do
! echo -n "abc_X($i,$j) = "
! for k in 1 2 3 4; do
! if [ $k == 4 ]; then
! echo "grid2(ii$i,ij$j,ik$k)*v$k"
! else
! echo -n "grid2(ii$i,ij$j,ik$k)*v$k + "
! fi
! done
! done
! done
! echo ""
! for j in 1 2 3 4; do
! echo -n "abc_X_Y($j) = "
! for k in 1 2 3 4; do
! if [ $k == 4 ]; then
! echo "abc_X($k,$j)*t$k"
! else
! echo -n "abc_X($k,$j)*t$k + "
! fi
! done
! done
! echo ""
! echo -n "val($ii) = "
! for k in 1 2 3 4; do
! if [ $k == 4 ]; then
! echo "abc_X_Y($k)*s$k"
! else
! echo -n "abc_X_Y($k)*s$k + "
! fi
! done
! echo ""
! done
t1 = t1d*dr1i
t2 = t2d*dr1i
t3 = t3d*dr1i
@ -1026,26 +995,29 @@ CONTAINS
v2 = v2o
v3 = v3o
v4 = v4o
abc_X(1, 1) = grid2(ii1, ij1, ik1)*v1 + grid2(ii1, ij1, ik2)*v2 + grid2(ii1, ij1, ik3)*v3 + grid2(ii1, ij1, ik4)*v4
abc_X(1, 2) = grid2(ii1, ij2, ik1)*v1 + grid2(ii1, ij2, ik2)*v2 + grid2(ii1, ij2, ik3)*v3 + grid2(ii1, ij2, ik4)*v4
abc_X(1, 3) = grid2(ii1, ij3, ik1)*v1 + grid2(ii1, ij3, ik2)*v2 + grid2(ii1, ij3, ik3)*v3 + grid2(ii1, ij3, ik4)*v4
abc_X(1, 4) = grid2(ii1, ij4, ik1)*v1 + grid2(ii1, ij4, ik2)*v2 + grid2(ii1, ij4, ik3)*v3 + grid2(ii1, ij4, ik4)*v4
abc_X(2, 1) = grid2(ii2, ij1, ik1)*v1 + grid2(ii2, ij1, ik2)*v2 + grid2(ii2, ij1, ik3)*v3 + grid2(ii2, ij1, ik4)*v4
abc_X(2, 2) = grid2(ii2, ij2, ik1)*v1 + grid2(ii2, ij2, ik2)*v2 + grid2(ii2, ij2, ik3)*v3 + grid2(ii2, ij2, ik4)*v4
abc_X(2, 3) = grid2(ii2, ij3, ik1)*v1 + grid2(ii2, ij3, ik2)*v2 + grid2(ii2, ij3, ik3)*v3 + grid2(ii2, ij3, ik4)*v4
abc_X(2, 4) = grid2(ii2, ij4, ik1)*v1 + grid2(ii2, ij4, ik2)*v2 + grid2(ii2, ij4, ik3)*v3 + grid2(ii2, ij4, ik4)*v4
abc_X(3, 1) = grid2(ii3, ij1, ik1)*v1 + grid2(ii3, ij1, ik2)*v2 + grid2(ii3, ij1, ik3)*v3 + grid2(ii3, ij1, ik4)*v4
abc_X(3, 2) = grid2(ii3, ij2, ik1)*v1 + grid2(ii3, ij2, ik2)*v2 + grid2(ii3, ij2, ik3)*v3 + grid2(ii3, ij2, ik4)*v4
abc_X(3, 3) = grid2(ii3, ij3, ik1)*v1 + grid2(ii3, ij3, ik2)*v2 + grid2(ii3, ij3, ik3)*v3 + grid2(ii3, ij3, ik4)*v4
abc_X(3, 4) = grid2(ii3, ij4, ik1)*v1 + grid2(ii3, ij4, ik2)*v2 + grid2(ii3, ij4, ik3)*v3 + grid2(ii3, ij4, ik4)*v4
abc_X(4, 1) = grid2(ii4, ij1, ik1)*v1 + grid2(ii4, ij1, ik2)*v2 + grid2(ii4, ij1, ik3)*v3 + grid2(ii4, ij1, ik4)*v4
abc_X(4, 2) = grid2(ii4, ij2, ik1)*v1 + grid2(ii4, ij2, ik2)*v2 + grid2(ii4, ij2, ik3)*v3 + grid2(ii4, ij2, ik4)*v4
abc_X(4, 3) = grid2(ii4, ij3, ik1)*v1 + grid2(ii4, ij3, ik2)*v2 + grid2(ii4, ij3, ik3)*v3 + grid2(ii4, ij3, ik4)*v4
abc_X(4, 4) = grid2(ii4, ij4, ik1)*v1 + grid2(ii4, ij4, ik2)*v2 + grid2(ii4, ij4, ik3)*v3 + grid2(ii4, ij4, ik4)*v4
abc_X(1, 1) = grid2(ii1, ij1, ik1)*v1 + grid2(ii1, ij1, ik2)*v2 + grid2(ii1, ij1, ik3)*v3 + grid2(ii1, ij1, ik4)*v4
abc_X(2, 1) = grid2(ii2, ij1, ik1)*v1 + grid2(ii2, ij1, ik2)*v2 + grid2(ii2, ij1, ik3)*v3 + grid2(ii2, ij1, ik4)*v4
abc_X(3, 1) = grid2(ii3, ij1, ik1)*v1 + grid2(ii3, ij1, ik2)*v2 + grid2(ii3, ij1, ik3)*v3 + grid2(ii3, ij1, ik4)*v4
abc_X(4, 1) = grid2(ii4, ij1, ik1)*v1 + grid2(ii4, ij1, ik2)*v2 + grid2(ii4, ij1, ik3)*v3 + grid2(ii4, ij1, ik4)*v4
abc_X_Y(1) = abc_X(1, 1)*t1 + abc_X(2, 1)*t2 + abc_X(3, 1)*t3 + abc_X(4, 1)*t4
abc_X(1, 2) = grid2(ii1, ij2, ik1)*v1 + grid2(ii1, ij2, ik2)*v2 + grid2(ii1, ij2, ik3)*v3 + grid2(ii1, ij2, ik4)*v4
abc_X(2, 2) = grid2(ii2, ij2, ik1)*v1 + grid2(ii2, ij2, ik2)*v2 + grid2(ii2, ij2, ik3)*v3 + grid2(ii2, ij2, ik4)*v4
abc_X(3, 2) = grid2(ii3, ij2, ik1)*v1 + grid2(ii3, ij2, ik2)*v2 + grid2(ii3, ij2, ik3)*v3 + grid2(ii3, ij2, ik4)*v4
abc_X(4, 2) = grid2(ii4, ij2, ik1)*v1 + grid2(ii4, ij2, ik2)*v2 + grid2(ii4, ij2, ik3)*v3 + grid2(ii4, ij2, ik4)*v4
abc_X_Y(2) = abc_X(1, 2)*t1 + abc_X(2, 2)*t2 + abc_X(3, 2)*t3 + abc_X(4, 2)*t4
abc_X(1, 3) = grid2(ii1, ij3, ik1)*v1 + grid2(ii1, ij3, ik2)*v2 + grid2(ii1, ij3, ik3)*v3 + grid2(ii1, ij3, ik4)*v4
abc_X(2, 3) = grid2(ii2, ij3, ik1)*v1 + grid2(ii2, ij3, ik2)*v2 + grid2(ii2, ij3, ik3)*v3 + grid2(ii2, ij3, ik4)*v4
abc_X(3, 3) = grid2(ii3, ij3, ik1)*v1 + grid2(ii3, ij3, ik2)*v2 + grid2(ii3, ij3, ik3)*v3 + grid2(ii3, ij3, ik4)*v4
abc_X(4, 3) = grid2(ii4, ij3, ik1)*v1 + grid2(ii4, ij3, ik2)*v2 + grid2(ii4, ij3, ik3)*v3 + grid2(ii4, ij3, ik4)*v4
abc_X_Y(3) = abc_X(1, 3)*t1 + abc_X(2, 3)*t2 + abc_X(3, 3)*t3 + abc_X(4, 3)*t4
abc_X(1, 4) = grid2(ii1, ij4, ik1)*v1 + grid2(ii1, ij4, ik2)*v2 + grid2(ii1, ij4, ik3)*v3 + grid2(ii1, ij4, ik4)*v4
abc_X(2, 4) = grid2(ii2, ij4, ik1)*v1 + grid2(ii2, ij4, ik2)*v2 + grid2(ii2, ij4, ik3)*v3 + grid2(ii2, ij4, ik4)*v4
abc_X(3, 4) = grid2(ii3, ij4, ik1)*v1 + grid2(ii3, ij4, ik2)*v2 + grid2(ii3, ij4, ik3)*v3 + grid2(ii3, ij4, ik4)*v4
abc_X(4, 4) = grid2(ii4, ij4, ik1)*v1 + grid2(ii4, ij4, ik2)*v2 + grid2(ii4, ij4, ik3)*v3 + grid2(ii4, ij4, ik4)*v4
abc_X_Y(4) = abc_X(1, 4)*t1 + abc_X(2, 4)*t2 + abc_X(3, 4)*t3 + abc_X(4, 4)*t4
val(1) = abc_X_Y(1)*s1 + abc_X_Y(2)*s2 + abc_X_Y(3)*s3 + abc_X_Y(4)*s4
@ -1058,26 +1030,6 @@ CONTAINS
s2 = s2d*dr2i
s3 = s3d*dr2i
s4 = s4d*dr2i
!! v1 = v1o
!! v2 = v2o
!! v3 = v3o
!! v4 = v4o
!! abc_X(1,1) = grid2(ii1,ij1,ik1)*v1 + grid2(ii1,ij1,ik2)*v2 + grid2(ii1,ij1,ik3)*v3 + grid2(ii1,ij1,ik4)*v4
!! abc_X(1,2) = grid2(ii1,ij2,ik1)*v1 + grid2(ii1,ij2,ik2)*v2 + grid2(ii1,ij2,ik3)*v3 + grid2(ii1,ij2,ik4)*v4
!! abc_X(1,3) = grid2(ii1,ij3,ik1)*v1 + grid2(ii1,ij3,ik2)*v2 + grid2(ii1,ij3,ik3)*v3 + grid2(ii1,ij3,ik4)*v4
!! abc_X(1,4) = grid2(ii1,ij4,ik1)*v1 + grid2(ii1,ij4,ik2)*v2 + grid2(ii1,ij4,ik3)*v3 + grid2(ii1,ij4,ik4)*v4
!! abc_X(2,1) = grid2(ii2,ij1,ik1)*v1 + grid2(ii2,ij1,ik2)*v2 + grid2(ii2,ij1,ik3)*v3 + grid2(ii2,ij1,ik4)*v4
!! abc_X(2,2) = grid2(ii2,ij2,ik1)*v1 + grid2(ii2,ij2,ik2)*v2 + grid2(ii2,ij2,ik3)*v3 + grid2(ii2,ij2,ik4)*v4
!! abc_X(2,3) = grid2(ii2,ij3,ik1)*v1 + grid2(ii2,ij3,ik2)*v2 + grid2(ii2,ij3,ik3)*v3 + grid2(ii2,ij3,ik4)*v4
!! abc_X(2,4) = grid2(ii2,ij4,ik1)*v1 + grid2(ii2,ij4,ik2)*v2 + grid2(ii2,ij4,ik3)*v3 + grid2(ii2,ij4,ik4)*v4
!! abc_X(3,1) = grid2(ii3,ij1,ik1)*v1 + grid2(ii3,ij1,ik2)*v2 + grid2(ii3,ij1,ik3)*v3 + grid2(ii3,ij1,ik4)*v4
!! abc_X(3,2) = grid2(ii3,ij2,ik1)*v1 + grid2(ii3,ij2,ik2)*v2 + grid2(ii3,ij2,ik3)*v3 + grid2(ii3,ij2,ik4)*v4
!! abc_X(3,3) = grid2(ii3,ij3,ik1)*v1 + grid2(ii3,ij3,ik2)*v2 + grid2(ii3,ij3,ik3)*v3 + grid2(ii3,ij3,ik4)*v4
!! abc_X(3,4) = grid2(ii3,ij4,ik1)*v1 + grid2(ii3,ij4,ik2)*v2 + grid2(ii3,ij4,ik3)*v3 + grid2(ii3,ij4,ik4)*v4
!! abc_X(4,1) = grid2(ii4,ij1,ik1)*v1 + grid2(ii4,ij1,ik2)*v2 + grid2(ii4,ij1,ik3)*v3 + grid2(ii4,ij1,ik4)*v4
!! abc_X(4,2) = grid2(ii4,ij2,ik1)*v1 + grid2(ii4,ij2,ik2)*v2 + grid2(ii4,ij2,ik3)*v3 + grid2(ii4,ij2,ik4)*v4
!! abc_X(4,3) = grid2(ii4,ij3,ik1)*v1 + grid2(ii4,ij3,ik2)*v2 + grid2(ii4,ij3,ik3)*v3 + grid2(ii4,ij3,ik4)*v4
!! abc_X(4,4) = grid2(ii4,ij4,ik1)*v1 + grid2(ii4,ij4,ik2)*v2 + grid2(ii4,ij4,ik3)*v3 + grid2(ii4,ij4,ik4)*v4
abc_X_Y(1) = abc_X(1, 1)*t1 + abc_X(2, 1)*t2 + abc_X(3, 1)*t3 + abc_X(4, 1)*t4
abc_X_Y(2) = abc_X(1, 2)*t1 + abc_X(2, 2)*t2 + abc_X(3, 2)*t3 + abc_X(4, 2)*t4
@ -1098,30 +1050,29 @@ CONTAINS
v2 = v2d*dr3i
v3 = v3d*dr3i
v4 = v4d*dr3i
abc_X(1, 1) = grid2(ii1, ij1, ik1)*v1 + grid2(ii1, ij1, ik2)*v2 + grid2(ii1, ij1, ik3)*v3 + grid2(ii1, ij1, ik4)*v4
abc_X(1, 2) = grid2(ii1, ij2, ik1)*v1 + grid2(ii1, ij2, ik2)*v2 + grid2(ii1, ij2, ik3)*v3 + grid2(ii1, ij2, ik4)*v4
abc_X(1, 3) = grid2(ii1, ij3, ik1)*v1 + grid2(ii1, ij3, ik2)*v2 + grid2(ii1, ij3, ik3)*v3 + grid2(ii1, ij3, ik4)*v4
abc_X(1, 4) = grid2(ii1, ij4, ik1)*v1 + grid2(ii1, ij4, ik2)*v2 + grid2(ii1, ij4, ik3)*v3 + grid2(ii1, ij4, ik4)*v4
abc_X(2, 1) = grid2(ii2, ij1, ik1)*v1 + grid2(ii2, ij1, ik2)*v2 + grid2(ii2, ij1, ik3)*v3 + grid2(ii2, ij1, ik4)*v4
abc_X(2, 2) = grid2(ii2, ij2, ik1)*v1 + grid2(ii2, ij2, ik2)*v2 + grid2(ii2, ij2, ik3)*v3 + grid2(ii2, ij2, ik4)*v4
abc_X(2, 3) = grid2(ii2, ij3, ik1)*v1 + grid2(ii2, ij3, ik2)*v2 + grid2(ii2, ij3, ik3)*v3 + grid2(ii2, ij3, ik4)*v4
abc_X(2, 4) = grid2(ii2, ij4, ik1)*v1 + grid2(ii2, ij4, ik2)*v2 + grid2(ii2, ij4, ik3)*v3 + grid2(ii2, ij4, ik4)*v4
abc_X(3, 1) = grid2(ii3, ij1, ik1)*v1 + grid2(ii3, ij1, ik2)*v2 + grid2(ii3, ij1, ik3)*v3 + grid2(ii3, ij1, ik4)*v4
abc_X(3, 2) = grid2(ii3, ij2, ik1)*v1 + grid2(ii3, ij2, ik2)*v2 + grid2(ii3, ij2, ik3)*v3 + grid2(ii3, ij2, ik4)*v4
abc_X(3, 3) = grid2(ii3, ij3, ik1)*v1 + grid2(ii3, ij3, ik2)*v2 + grid2(ii3, ij3, ik3)*v3 + grid2(ii3, ij3, ik4)*v4
abc_X(3, 4) = grid2(ii3, ij4, ik1)*v1 + grid2(ii3, ij4, ik2)*v2 + grid2(ii3, ij4, ik3)*v3 + grid2(ii3, ij4, ik4)*v4
abc_X(4, 1) = grid2(ii4, ij1, ik1)*v1 + grid2(ii4, ij1, ik2)*v2 + grid2(ii4, ij1, ik3)*v3 + grid2(ii4, ij1, ik4)*v4
abc_X(4, 2) = grid2(ii4, ij2, ik1)*v1 + grid2(ii4, ij2, ik2)*v2 + grid2(ii4, ij2, ik3)*v3 + grid2(ii4, ij2, ik4)*v4
abc_X(4, 3) = grid2(ii4, ij3, ik1)*v1 + grid2(ii4, ij3, ik2)*v2 + grid2(ii4, ij3, ik3)*v3 + grid2(ii4, ij3, ik4)*v4
abc_X(4, 4) = grid2(ii4, ij4, ik1)*v1 + grid2(ii4, ij4, ik2)*v2 + grid2(ii4, ij4, ik3)*v3 + grid2(ii4, ij4, ik4)*v4
abc_X(1, 1) = grid2(ii1, ij1, ik1)*v1 + grid2(ii1, ij1, ik2)*v2 + grid2(ii1, ij1, ik3)*v3 + grid2(ii1, ij1, ik4)*v4
abc_X(2, 1) = grid2(ii2, ij1, ik1)*v1 + grid2(ii2, ij1, ik2)*v2 + grid2(ii2, ij1, ik3)*v3 + grid2(ii2, ij1, ik4)*v4
abc_X(3, 1) = grid2(ii3, ij1, ik1)*v1 + grid2(ii3, ij1, ik2)*v2 + grid2(ii3, ij1, ik3)*v3 + grid2(ii3, ij1, ik4)*v4
abc_X(4, 1) = grid2(ii4, ij1, ik1)*v1 + grid2(ii4, ij1, ik2)*v2 + grid2(ii4, ij1, ik3)*v3 + grid2(ii4, ij1, ik4)*v4
abc_X_Y(1) = abc_X(1, 1)*t1 + abc_X(2, 1)*t2 + abc_X(3, 1)*t3 + abc_X(4, 1)*t4
abc_X(1, 2) = grid2(ii1, ij2, ik1)*v1 + grid2(ii1, ij2, ik2)*v2 + grid2(ii1, ij2, ik3)*v3 + grid2(ii1, ij2, ik4)*v4
abc_X(2, 2) = grid2(ii2, ij2, ik1)*v1 + grid2(ii2, ij2, ik2)*v2 + grid2(ii2, ij2, ik3)*v3 + grid2(ii2, ij2, ik4)*v4
abc_X(3, 2) = grid2(ii3, ij2, ik1)*v1 + grid2(ii3, ij2, ik2)*v2 + grid2(ii3, ij2, ik3)*v3 + grid2(ii3, ij2, ik4)*v4
abc_X(4, 2) = grid2(ii4, ij2, ik1)*v1 + grid2(ii4, ij2, ik2)*v2 + grid2(ii4, ij2, ik3)*v3 + grid2(ii4, ij2, ik4)*v4
abc_X_Y(2) = abc_X(1, 2)*t1 + abc_X(2, 2)*t2 + abc_X(3, 2)*t3 + abc_X(4, 2)*t4
abc_X(1, 3) = grid2(ii1, ij3, ik1)*v1 + grid2(ii1, ij3, ik2)*v2 + grid2(ii1, ij3, ik3)*v3 + grid2(ii1, ij3, ik4)*v4
abc_X(2, 3) = grid2(ii2, ij3, ik1)*v1 + grid2(ii2, ij3, ik2)*v2 + grid2(ii2, ij3, ik3)*v3 + grid2(ii2, ij3, ik4)*v4
abc_X(3, 3) = grid2(ii3, ij3, ik1)*v1 + grid2(ii3, ij3, ik2)*v2 + grid2(ii3, ij3, ik3)*v3 + grid2(ii3, ij3, ik4)*v4
abc_X(4, 3) = grid2(ii4, ij3, ik1)*v1 + grid2(ii4, ij3, ik2)*v2 + grid2(ii4, ij3, ik3)*v3 + grid2(ii4, ij3, ik4)*v4
abc_X_Y(3) = abc_X(1, 3)*t1 + abc_X(2, 3)*t2 + abc_X(3, 3)*t3 + abc_X(4, 3)*t4
abc_X(1, 4) = grid2(ii1, ij4, ik1)*v1 + grid2(ii1, ij4, ik2)*v2 + grid2(ii1, ij4, ik3)*v3 + grid2(ii1, ij4, ik4)*v4
abc_X(2, 4) = grid2(ii2, ij4, ik1)*v1 + grid2(ii2, ij4, ik2)*v2 + grid2(ii2, ij4, ik3)*v3 + grid2(ii2, ij4, ik4)*v4
abc_X(3, 4) = grid2(ii3, ij4, ik1)*v1 + grid2(ii3, ij4, ik2)*v2 + grid2(ii3, ij4, ik3)*v3 + grid2(ii3, ij4, ik4)*v4
abc_X(4, 4) = grid2(ii4, ij4, ik1)*v1 + grid2(ii4, ij4, ik2)*v2 + grid2(ii4, ij4, ik3)*v3 + grid2(ii4, ij4, ik4)*v4
abc_X_Y(4) = abc_X(1, 4)*t1 + abc_X(2, 4)*t2 + abc_X(3, 4)*t3 + abc_X(4, 4)*t4
val(3) = abc_X_Y(1)*s1 + abc_X_Y(2)*s2 + abc_X_Y(3)*s3 + abc_X_Y(4)*s4
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
fac = grid(i, j, k)
ft1 = ft1 + val(1)*fac
@ -1141,6 +1092,7 @@ CONTAINS
Forces(2, LIndMM) = Forces(2, LIndMM) + LForces(2, LindMM)
Forces(3, LIndMM) = Forces(3, LIndMM) + LForces(3, LindMM)
END DO Atoms
!$OMP END PARALLEL DO
END DO Radius
!
! Debug Statement
@ -1214,7 +1166,7 @@ CONTAINS
LOGICAL :: shells
CHARACTER(len=*), PARAMETER :: routineN = 'qmmm_forces_with_gaussian_LR', &
routineP = moduleN//':'//routineN
routineP = moduleN//':'//routineN
INTEGER :: handle, i, Imm, IndMM, IRadTyp, ix, j, &
k, LIndMM, my_i, my_j, my_k, myind, &
@ -1248,9 +1200,16 @@ CONTAINS
pot => potentials(IRadTyp)%pot
dx = Pot%dx
pot0_2 => Pot%pot0_2
!$OMP PARALLEL DO DEFAULT(NONE) &
!$OMP SHARED(pot, par_scheme, para_env, dvol, mm_atom_index, mm_particles, dOmmOqm) &
!$OMP SHARED(mm_cell, mm_charges, dx, LForces, Forces, qmmm_spherical_cutoff, shells, dr1, dr2, dr3, gbo, bo) &
!$OMP SHARED(IRadTyp, pot0_2, grid) &
!$OMP PRIVATE(Imm, myind, ra, LIndMM, IndMM, qt, rt1, rt2, rt3, ft1, ft2, ft3, i, j, k, sph_chrg_factor) &
!$OMP PRIVATE(my_k, my_j, my_i, xs3, xs2, xs1, rv1, rv2, rv3, r, ix, rx, rx2, r2, Term, fac) &
!$OMP PRIVATE(rd1, rd2, rd3)
Atoms: DO Imm = 1, SIZE(pot%mm_atom_index)
IF (par_scheme == do_par_atom) THEN
myind = myind + 1
myind = Imm + (IRadTyp - 1)*SIZE(pot%mm_atom_index)
IF (MOD(myind, para_env%num_pe) /= para_env%mepos) CYCLE
END IF
LIndMM = pot%mm_atom_index(Imm)
@ -1319,6 +1278,7 @@ CONTAINS
Forces(2, LIndMM) = Forces(2, LIndMM) + LForces(2, LindMM)
Forces(3, LIndMM) = Forces(3, LIndMM) + LForces(3, LindMM)
END DO Atoms
!$OMP END PARALLEL DO
END DO Radius
!
! Debug Statement
@ -1378,7 +1338,7 @@ CONTAINS
TYPE(cell_type), POINTER :: mm_cell
CHARACTER(len=*), PARAMETER :: routineN = 'qmmm_debug_forces', &
routineP = moduleN//':'//routineN
routineP = moduleN//':'//routineN
INTEGER :: handle, I, IndMM, iw, J, K
REAL(KIND=dp) :: Coord_save
@ -1526,7 +1486,7 @@ CONTAINS
INTEGER, INTENT(IN) :: iw
CHARACTER(len=*), PARAMETER :: routineN = 'debug_integrate_gf_rspace_NoPBC', &
routineP = moduleN//':'//routineN
routineP = moduleN//':'//routineN
INTEGER :: handle, i, igrid, k, ngrids
INTEGER, DIMENSION(2, 3) :: bo2
@ -1648,7 +1608,7 @@ CONTAINS
LOGICAL :: shells
CHARACTER(len=*), PARAMETER :: routineN = 'debug_qmmm_forces_with_gauss_LG', &
routineP = moduleN//':'//routineN
routineP = moduleN//':'//routineN
INTEGER :: handle, I, igrid, IndMM, J, K, ngrids
REAL(KIND=dp) :: Coord_save
@ -1764,7 +1724,7 @@ CONTAINS
LOGICAL :: shells
CHARACTER(len=*), PARAMETER :: routineN = 'debug_qmmm_forces_with_gauss_LR', &
routineP = moduleN//':'//routineN
routineP = moduleN//':'//routineN
INTEGER :: handle, I, igrid, IndMM, J, K, ngrids
REAL(KIND=dp) :: Coord_save