diff --git a/src/qmmm_gpw_energy.F b/src/qmmm_gpw_energy.F index feb9e5d2a9..24433746aa 100644 --- a/src/qmmm_gpw_energy.F +++ b/src/qmmm_gpw_energy.F @@ -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 diff --git a/src/qmmm_gpw_forces.F b/src/qmmm_gpw_forces.F index 64760f446b..8abe0c60af 100644 --- a/src/qmmm_gpw_forces.F +++ b/src/qmmm_gpw_forces.F @@ -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