1!--------------------------------------------------------------------------------------------------! 2! CP2K: A general program to perform molecular dynamics simulations ! 3! Copyright (C) 2000 - 2020 CP2K developers group ! 4!--------------------------------------------------------------------------------------------------! 5 6! ************************************************************************************************** 7!> \par History 8!> Efficient tersoff implementation and general "lifting" of manybody_potential module 9!> 12.2007 [tlaino] - Splitting manybody module : In this module we should only 10!> keep the main routines for computing energy and forces of 11!> manybody potentials. Each potential should have his own module! 12!> \author CJM, I-Feng W. Kuo, Teodoro Laino 13! ************************************************************************************************** 14MODULE manybody_potential 15 16 USE atomic_kind_types, ONLY: atomic_kind_type 17 USE cell_types, ONLY: cell_type 18 USE cp_para_types, ONLY: cp_para_env_type 19 USE distribution_1d_types, ONLY: distribution_1d_type 20 USE fist_neighbor_list_types, ONLY: fist_neighbor_type,& 21 neighbor_kind_pairs_type 22 USE fist_nonbond_env_types, ONLY: eam_type,& 23 fist_nonbond_env_get,& 24 fist_nonbond_env_type,& 25 pos_type 26 USE input_section_types, ONLY: section_vals_type 27 USE kinds, ONLY: dp 28 USE manybody_eam, ONLY: get_force_eam 29 USE manybody_quip, ONLY: quip_add_force_virial,& 30 quip_energy_store_force_virial 31 USE manybody_siepmann, ONLY: destroy_siepmann_arrays,& 32 print_nr_ions_siepmann,& 33 setup_siepmann_arrays,& 34 siepmann_energy,& 35 siepmann_forces_v2,& 36 siepmann_forces_v3 37 USE manybody_tersoff, ONLY: destroy_tersoff_arrays,& 38 setup_tersoff_arrays,& 39 tersoff_energy,& 40 tersoff_forces 41 USE mathlib, ONLY: matvec_3x3 42 USE message_passing, ONLY: mp_sum 43 USE pair_potential_types, ONLY: & 44 ea_type, eam_pot_type, pair_potential_pp_type, pair_potential_single_type, quip_type, & 45 siepmann_pot_type, siepmann_type, tersoff_pot_type, tersoff_type 46 USE particle_types, ONLY: particle_type 47 USE util, ONLY: sort 48#include "./base/base_uses.f90" 49 50 IMPLICIT NONE 51 52 PRIVATE 53 PUBLIC :: energy_manybody 54 PUBLIC :: force_nonbond_manybody 55 CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'manybody_potential' 56 57CONTAINS 58 59! ************************************************************************************************** 60!> \brief computes the embedding contribution to the energy 61!> \param fist_nonbond_env ... 62!> \param atomic_kind_set ... 63!> \param local_particles ... 64!> \param particle_set ... 65!> \param cell ... 66!> \param pot_manybody ... 67!> \param para_env ... 68!> \param mm_section ... 69!> \par History 70!> tlaino [2007] - New algorithm for tersoff potential 71!> \author CJM, I-Feng W. Kuo, Teodoro Laino 72! ************************************************************************************************** 73 SUBROUTINE energy_manybody(fist_nonbond_env, atomic_kind_set, local_particles, & 74 particle_set, cell, pot_manybody, para_env, mm_section) 75 76 TYPE(fist_nonbond_env_type), POINTER :: fist_nonbond_env 77 TYPE(atomic_kind_type), POINTER :: atomic_kind_set(:) 78 TYPE(distribution_1d_type), POINTER :: local_particles 79 TYPE(particle_type), POINTER :: particle_set(:) 80 TYPE(cell_type), POINTER :: cell 81 REAL(dp), INTENT(INOUT) :: pot_manybody 82 TYPE(cp_para_env_type), OPTIONAL, POINTER :: para_env 83 TYPE(section_vals_type), POINTER :: mm_section 84 85 CHARACTER(LEN=*), PARAMETER :: routineN = 'energy_manybody', & 86 routineP = moduleN//':'//routineN 87 88 INTEGER :: atom_a, atom_b, handle, i, iend, ifirst, igrp, ikind, ilast, ilist, indexa, & 89 ipair, iparticle, iparticle_local, istart, iunique, jkind, junique, mpair, nkinds, & 90 nloc_size, npairs, nparticle, nparticle_local, nr_h3O, nr_o, nr_oh, nunique 91 INTEGER, DIMENSION(:), POINTER :: glob_loc_list_a, work_list 92 INTEGER, DIMENSION(:, :), POINTER :: glob_loc_list, list, sort_list 93 LOGICAL :: any_quip, any_siepmann, any_tersoff 94 REAL(KIND=dp) :: drij, embed, pot_loc, pot_quip, qr, & 95 rab2_max, rij(3) 96 REAL(KIND=dp), DIMENSION(3) :: cell_v, cvi 97 REAL(KIND=dp), DIMENSION(:, :), POINTER :: glob_cell_v 98 REAL(KIND=dp), POINTER :: fembed(:) 99 TYPE(eam_pot_type), POINTER :: eam 100 TYPE(eam_type), DIMENSION(:), POINTER :: eam_data 101 TYPE(fist_neighbor_type), POINTER :: nonbonded 102 TYPE(neighbor_kind_pairs_type), POINTER :: neighbor_kind_pair 103 TYPE(pair_potential_pp_type), POINTER :: potparm 104 TYPE(pair_potential_single_type), POINTER :: pot 105 TYPE(pos_type), DIMENSION(:), POINTER :: r_last_update_pbc 106 TYPE(siepmann_pot_type), POINTER :: siepmann 107 TYPE(tersoff_pot_type), POINTER :: tersoff 108 109 NULLIFY (eam, siepmann, tersoff) 110 any_tersoff = .FALSE. 111 any_siepmann = .FALSE. 112 any_quip = .FALSE. 113 CALL timeset(routineN, handle) 114 CALL fist_nonbond_env_get(fist_nonbond_env, r_last_update_pbc=r_last_update_pbc, & 115 potparm=potparm, eam_data=eam_data) 116 ! EAM requires a single loop 117 DO ikind = 1, SIZE(atomic_kind_set) 118 pot => potparm%pot(ikind, ikind)%pot 119 DO i = 1, SIZE(pot%type) 120 IF (pot%type(i) /= ea_type) CYCLE 121 eam => pot%set(i)%eam 122 nparticle = SIZE(particle_set) 123 ALLOCATE (fembed(nparticle)) 124 fembed(:) = 0._dp 125 CPASSERT(ASSOCIATED(eam_data)) 126 ! computation of embedding function and energy 127 nparticle_local = local_particles%n_el(ikind) 128 DO iparticle_local = 1, nparticle_local 129 iparticle = local_particles%list(ikind)%array(iparticle_local) 130 indexa = INT(eam_data(iparticle)%rho/eam%drhoar) + 1 131 IF (indexa > eam%npoints - 1) indexa = eam%npoints - 1 132 qr = eam_data(iparticle)%rho - eam%rhoval(indexa) 133 134 embed = eam%frho(indexa) + qr*eam%frhop(indexa) 135 fembed(iparticle) = eam%frhop(indexa) + qr*(eam%frhop(indexa + 1) - eam%frhop(indexa))/eam%drhoar 136 137 pot_manybody = pot_manybody + embed 138 END DO 139 ! communicate data 140 CALL mp_sum(fembed, para_env%group) 141 DO iparticle = 1, nparticle 142 IF (particle_set(iparticle)%atomic_kind%kind_number == ikind) THEN 143 eam_data(iparticle)%f_embed = fembed(iparticle) 144 END IF 145 END DO 146 147 DEALLOCATE (fembed) 148 END DO 149 END DO 150 ! Other manybody potential 151 DO ikind = 1, SIZE(atomic_kind_set) 152 DO jkind = ikind, SIZE(atomic_kind_set) 153 pot => potparm%pot(ikind, jkind)%pot 154 any_tersoff = any_tersoff .OR. ANY(pot%type == tersoff_type) 155 any_quip = any_quip .OR. ANY(pot%type == quip_type) 156 any_siepmann = any_siepmann .OR. ANY(pot%type == siepmann_type) 157 END DO 158 END DO 159 CALL fist_nonbond_env_get(fist_nonbond_env, nonbonded=nonbonded, natom_types=nkinds) 160 ! QUIP 161 IF (any_quip) THEN 162 CALL quip_energy_store_force_virial(particle_set, cell, atomic_kind_set, potparm, & 163 fist_nonbond_env, pot_quip, para_env) 164 pot_manybody = pot_manybody + pot_quip 165 ENDIF 166 ! TERSOFF 167 IF (any_tersoff) THEN 168 NULLIFY (glob_loc_list, glob_cell_v, glob_loc_list_a) 169 CALL setup_tersoff_arrays(nonbonded, potparm, glob_loc_list, glob_cell_v, glob_loc_list_a, cell) 170 DO ilist = 1, nonbonded%nlists 171 neighbor_kind_pair => nonbonded%neighbor_kind_pairs(ilist) 172 npairs = neighbor_kind_pair%npairs 173 IF (npairs == 0) CYCLE 174 Kind_Group_Loop: DO igrp = 1, neighbor_kind_pair%ngrp_kind 175 istart = neighbor_kind_pair%grp_kind_start(igrp) 176 iend = neighbor_kind_pair%grp_kind_end(igrp) 177 ikind = neighbor_kind_pair%ij_kind(1, igrp) 178 jkind = neighbor_kind_pair%ij_kind(2, igrp) 179 list => neighbor_kind_pair%list 180 cvi = neighbor_kind_pair%cell_vector 181 pot => potparm%pot(ikind, jkind)%pot 182 DO i = 1, SIZE(pot%type) 183 IF (pot%type(i) /= tersoff_type) CYCLE 184 rab2_max = pot%set(i)%tersoff%rcutsq 185 CALL matvec_3x3(cell_v, cell%hmat, cvi) 186 pot => potparm%pot(ikind, jkind)%pot 187 tersoff => pot%set(i)%tersoff 188 npairs = iend - istart + 1 189 IF (npairs /= 0) THEN 190 ALLOCATE (sort_list(2, npairs), work_list(npairs)) 191 sort_list = list(:, istart:iend) 192 ! Sort the list of neighbors, this increases the efficiency for single 193 ! potential contributions 194 CALL sort(sort_list(1, :), npairs, work_list) 195 DO ipair = 1, npairs 196 work_list(ipair) = sort_list(2, work_list(ipair)) 197 END DO 198 sort_list(2, :) = work_list 199 ! find number of unique elements of array index 1 200 nunique = 1 201 DO ipair = 1, npairs - 1 202 IF (sort_list(1, ipair + 1) /= sort_list(1, ipair)) nunique = nunique + 1 203 END DO 204 ipair = 1 205 junique = sort_list(1, ipair) 206 ifirst = 1 207 DO iunique = 1, nunique 208 atom_a = junique 209 IF (glob_loc_list_a(ifirst) > atom_a) CYCLE 210 DO mpair = ifirst, SIZE(glob_loc_list_a) 211 IF (glob_loc_list_a(mpair) == atom_a) EXIT 212 END DO 213 ifirst = mpair 214 DO mpair = ifirst, SIZE(glob_loc_list_a) 215 IF (glob_loc_list_a(mpair) /= atom_a) EXIT 216 END DO 217 ilast = mpair - 1 218 nloc_size = 0 219 IF (ifirst /= 0) nloc_size = ilast - ifirst + 1 220 DO WHILE (ipair <= npairs) 221 IF (sort_list(1, ipair) /= junique) EXIT 222 atom_b = sort_list(2, ipair) 223 ! Energy terms 224 pot_loc = 0.0_dp 225 rij(:) = r_last_update_pbc(atom_b)%r(:) - r_last_update_pbc(atom_a)%r(:) + cell_v 226 drij = DOT_PRODUCT(rij, rij) 227 ipair = ipair + 1 228 IF (drij > rab2_max) CYCLE 229 drij = SQRT(drij) 230 CALL tersoff_energy(pot_loc, tersoff, r_last_update_pbc, atom_a, atom_b, nloc_size, & 231 glob_loc_list(:, ifirst:ilast), glob_cell_v(:, ifirst:ilast), cell_v, drij) 232 pot_manybody = pot_manybody + 0.5_dp*pot_loc 233 END DO 234 ifirst = ilast + 1 235 IF (ipair <= npairs) junique = sort_list(1, ipair) 236 END DO 237 DEALLOCATE (sort_list, work_list) 238 END IF 239 END DO 240 END DO Kind_Group_Loop 241 END DO 242 CALL destroy_tersoff_arrays(glob_loc_list, glob_cell_v, glob_loc_list_a) 243 END IF 244 245 !SIEPMANN POTENTIAL 246 IF (any_siepmann) THEN 247 NULLIFY (glob_loc_list, glob_cell_v, glob_loc_list_a) 248 nr_oh = 0 249 nr_h3O = 0 250 nr_o = 0 251 CALL setup_siepmann_arrays(nonbonded, potparm, glob_loc_list, glob_cell_v, glob_loc_list_a, cell) 252 DO ilist = 1, nonbonded%nlists 253 neighbor_kind_pair => nonbonded%neighbor_kind_pairs(ilist) 254 npairs = neighbor_kind_pair%npairs 255 IF (npairs == 0) CYCLE 256 Kind_Group_Loop_2: DO igrp = 1, neighbor_kind_pair%ngrp_kind 257 istart = neighbor_kind_pair%grp_kind_start(igrp) 258 iend = neighbor_kind_pair%grp_kind_end(igrp) 259 ikind = neighbor_kind_pair%ij_kind(1, igrp) 260 jkind = neighbor_kind_pair%ij_kind(2, igrp) 261 list => neighbor_kind_pair%list 262 cvi = neighbor_kind_pair%cell_vector 263 pot => potparm%pot(ikind, jkind)%pot 264 DO i = 1, SIZE(pot%type) 265 IF (pot%type(i) /= siepmann_type) CYCLE 266 rab2_max = pot%set(i)%siepmann%rcutsq 267 CALL matvec_3x3(cell_v, cell%hmat, cvi) 268 pot => potparm%pot(ikind, jkind)%pot 269 siepmann => pot%set(i)%siepmann 270 npairs = iend - istart + 1 271 IF (npairs /= 0) THEN 272 ALLOCATE (sort_list(2, npairs), work_list(npairs)) 273 sort_list = list(:, istart:iend) 274 ! Sort the list of neighbors, this increases the efficiency for single 275 ! potential contributions 276 CALL sort(sort_list(1, :), npairs, work_list) 277 DO ipair = 1, npairs 278 work_list(ipair) = sort_list(2, work_list(ipair)) 279 END DO 280 sort_list(2, :) = work_list 281 ! find number of unique elements of array index 1 282 nunique = 1 283 DO ipair = 1, npairs - 1 284 IF (sort_list(1, ipair + 1) /= sort_list(1, ipair)) nunique = nunique + 1 285 END DO 286 ipair = 1 287 junique = sort_list(1, ipair) 288 ifirst = 1 289 DO iunique = 1, nunique 290 atom_a = junique 291 IF (glob_loc_list_a(ifirst) > atom_a) CYCLE 292 DO mpair = ifirst, SIZE(glob_loc_list_a) 293 IF (glob_loc_list_a(mpair) == atom_a) EXIT 294 END DO 295 ifirst = mpair 296 DO mpair = ifirst, SIZE(glob_loc_list_a) 297 IF (glob_loc_list_a(mpair) /= atom_a) EXIT 298 END DO 299 ilast = mpair - 1 300 nloc_size = 0 301 IF (ifirst /= 0) nloc_size = ilast - ifirst + 1 302 DO WHILE (ipair <= npairs) 303 IF (sort_list(1, ipair) /= junique) EXIT 304 atom_b = sort_list(2, ipair) 305 ! Energy terms 306 pot_loc = 0.0_dp 307 rij(:) = r_last_update_pbc(atom_b)%r(:) - r_last_update_pbc(atom_a)%r(:) + cell_v 308 drij = DOT_PRODUCT(rij, rij) 309 ipair = ipair + 1 310 IF (drij > rab2_max) CYCLE 311 drij = SQRT(drij) 312 CALL siepmann_energy(pot_loc, siepmann, r_last_update_pbc, atom_a, atom_b, nloc_size, & 313 glob_loc_list(:, ifirst:ilast), cell_v, cell, drij, & 314 particle_set, nr_oh, nr_h3O, nr_o) 315 pot_manybody = pot_manybody + pot_loc 316 END DO 317 ifirst = ilast + 1 318 IF (ipair <= npairs) junique = sort_list(1, ipair) 319 END DO 320 DEALLOCATE (sort_list, work_list) 321 END IF 322 END DO 323 END DO Kind_Group_Loop_2 324 END DO 325 CALL destroy_siepmann_arrays(glob_loc_list, glob_cell_v, glob_loc_list_a) 326 CALL print_nr_ions_siepmann(nr_oh, mm_section, para_env, print_oh=.TRUE., & 327 print_h3o=.FALSE., print_o=.FALSE.) 328 CALL print_nr_ions_siepmann(nr_h3o, mm_section, para_env, print_oh=.FALSE., & 329 print_h3o=.TRUE., print_o=.FALSE.) 330 CALL print_nr_ions_siepmann(nr_o, mm_section, para_env, print_oh=.FALSE., & 331 print_h3o=.FALSE., print_o=.TRUE.) 332 END IF 333 334 CALL timestop(handle) 335 END SUBROUTINE energy_manybody 336 337! ************************************************************************************************** 338!> \brief ... 339!> \param fist_nonbond_env ... 340!> \param particle_set ... 341!> \param cell ... 342!> \param f_nonbond ... 343!> \param pv_nonbond ... 344!> \param use_virial ... 345!> \par History 346!> Fast implementation of the tersoff potential - [tlaino] 2007 347!> \author I-Feng W. Kuo, Teodoro Laino 348! ************************************************************************************************** 349 SUBROUTINE force_nonbond_manybody(fist_nonbond_env, particle_set, cell, & 350 f_nonbond, pv_nonbond, use_virial) 351 352 TYPE(fist_nonbond_env_type), POINTER :: fist_nonbond_env 353 TYPE(particle_type), DIMENSION(:), POINTER :: particle_set 354 TYPE(cell_type), POINTER :: cell 355 REAL(KIND=dp), DIMENSION(:, :), INTENT(INOUT) :: f_nonbond, pv_nonbond 356 LOGICAL, INTENT(IN) :: use_virial 357 358 CHARACTER(LEN=*), PARAMETER :: routineN = 'force_nonbond_manybody', & 359 routineP = moduleN//':'//routineN 360 361 INTEGER :: atom_a, atom_b, handle, i, i_a, i_b, iend, ifirst, igrp, ikind, ilast, ilist, & 362 ipair, istart, iunique, jkind, junique, kind_a, kind_b, mpair, nkinds, nloc_size, npairs, & 363 nunique 364 INTEGER, ALLOCATABLE, DIMENSION(:, :) :: eam_kinds_index 365 INTEGER, DIMENSION(:), POINTER :: glob_loc_list_a, work_list 366 INTEGER, DIMENSION(:, :), POINTER :: glob_loc_list, list, sort_list 367 LOGICAL :: any_quip, any_siepmann, any_tersoff 368 REAL(KIND=dp) :: f_eam, fac, fr(3), ptens11, ptens12, ptens13, ptens21, ptens22, ptens23, & 369 ptens31, ptens32, ptens33, rab(3), rab2, rab2_max, rtmp(3) 370 REAL(KIND=dp), DIMENSION(3) :: cell_v, cvi 371 REAL(KIND=dp), DIMENSION(:, :), POINTER :: glob_cell_v 372 TYPE(eam_pot_type), POINTER :: eam_a, eam_b 373 TYPE(eam_type), DIMENSION(:), POINTER :: eam_data 374 TYPE(fist_neighbor_type), POINTER :: nonbonded 375 TYPE(neighbor_kind_pairs_type), POINTER :: neighbor_kind_pair 376 TYPE(pair_potential_pp_type), POINTER :: potparm 377 TYPE(pair_potential_single_type), POINTER :: pot 378 TYPE(pos_type), DIMENSION(:), POINTER :: r_last_update_pbc 379 TYPE(siepmann_pot_type), POINTER :: siepmann 380 TYPE(tersoff_pot_type), POINTER :: tersoff 381 382 any_tersoff = .FALSE. 383 any_quip = .FALSE. 384 any_siepmann = .FALSE. 385 CALL timeset(routineN, handle) 386 NULLIFY (eam_a, eam_b, tersoff, siepmann) 387 388 CALL fist_nonbond_env_get(fist_nonbond_env, nonbonded=nonbonded, potparm=potparm, & 389 natom_types=nkinds, eam_data=eam_data, r_last_update_pbc=r_last_update_pbc) 390 391 ! Initializing the potential energy, pressure tensor and force 392 IF (use_virial) THEN 393 ptens11 = 0.0_dp; ptens12 = 0.0_dp; ptens13 = 0.0_dp 394 ptens21 = 0.0_dp; ptens22 = 0.0_dp; ptens23 = 0.0_dp 395 ptens31 = 0.0_dp; ptens32 = 0.0_dp; ptens33 = 0.0_dp 396 END IF 397 398 nkinds = SIZE(potparm%pot, 1) 399 ALLOCATE (eam_kinds_index(nkinds, nkinds)) 400 eam_kinds_index = -1 401 DO ikind = 1, nkinds 402 DO jkind = ikind, nkinds 403 DO i = 1, SIZE(potparm%pot(ikind, jkind)%pot%type) 404 IF (potparm%pot(ikind, jkind)%pot%type(i) == ea_type) THEN 405 ! At the moment we allow only 1 EAM per each kinds pair.. 406 CPASSERT(eam_kinds_index(ikind, jkind) == -1) 407 CPASSERT(eam_kinds_index(jkind, ikind) == -1) 408 eam_kinds_index(ikind, jkind) = i 409 eam_kinds_index(jkind, ikind) = i 410 END IF 411 END DO 412 END DO 413 END DO 414 415 ! QUIP 416 IF (use_virial) & 417 CALL quip_add_force_virial(fist_nonbond_env, f_nonbond, pv_nonbond) 418 419 ! starting the force loop 420 DO ilist = 1, nonbonded%nlists 421 neighbor_kind_pair => nonbonded%neighbor_kind_pairs(ilist) 422 npairs = neighbor_kind_pair%npairs 423 IF (npairs == 0) CYCLE 424 Kind_Group_Loop1: DO igrp = 1, neighbor_kind_pair%ngrp_kind 425 istart = neighbor_kind_pair%grp_kind_start(igrp) 426 iend = neighbor_kind_pair%grp_kind_end(igrp) 427 ikind = neighbor_kind_pair%ij_kind(1, igrp) 428 jkind = neighbor_kind_pair%ij_kind(2, igrp) 429 list => neighbor_kind_pair%list 430 cvi = neighbor_kind_pair%cell_vector 431 pot => potparm%pot(ikind, jkind)%pot 432 IF (pot%no_mb) CYCLE 433 rab2_max = pot%rcutsq 434 CALL matvec_3x3(cell_v, cell%hmat, cvi) 435 any_tersoff = any_tersoff .OR. ANY(pot%type == tersoff_type) 436 any_siepmann = any_siepmann .OR. ANY(pot%type == siepmann_type) 437 i = eam_kinds_index(ikind, jkind) 438 IF (i == -1) CYCLE 439 ! EAM 440 CPASSERT(ASSOCIATED(eam_data)) 441 DO ipair = istart, iend 442 atom_a = list(1, ipair) 443 atom_b = list(2, ipair) 444 fac = 1.0_dp 445 IF (atom_a == atom_b) fac = 0.5_dp 446 kind_a = particle_set(atom_a)%atomic_kind%kind_number 447 kind_b = particle_set(atom_b)%atomic_kind%kind_number 448 i_a = eam_kinds_index(kind_a, kind_a) 449 i_b = eam_kinds_index(kind_b, kind_b) 450 eam_a => potparm%pot(kind_a, kind_a)%pot%set(i_a)%eam 451 eam_b => potparm%pot(kind_b, kind_b)%pot%set(i_b)%eam 452 453 !set this outside the potential type in case need multiple potentials 454 !Do everything necessary for EAM here 455 rab = r_last_update_pbc(atom_b)%r - r_last_update_pbc(atom_a)%r 456 rab = rab + cell_v 457 rab2 = rab(1)*rab(1) + rab(2)*rab(2) + rab(3)*rab(3) 458 IF (rab2 <= rab2_max) THEN 459 CALL get_force_eam(rab2, eam_a, eam_b, eam_data, atom_a, atom_b, f_eam) 460 f_eam = f_eam*fac 461 462 fr(1) = -f_eam*rab(1) 463 fr(2) = -f_eam*rab(2) 464 fr(3) = -f_eam*rab(3) 465 f_nonbond(1, atom_a) = f_nonbond(1, atom_a) - fr(1) 466 f_nonbond(2, atom_a) = f_nonbond(2, atom_a) - fr(2) 467 f_nonbond(3, atom_a) = f_nonbond(3, atom_a) - fr(3) 468 469 f_nonbond(1, atom_b) = f_nonbond(1, atom_b) + fr(1) 470 f_nonbond(2, atom_b) = f_nonbond(2, atom_b) + fr(2) 471 f_nonbond(3, atom_b) = f_nonbond(3, atom_b) + fr(3) 472 IF (use_virial) THEN 473 ptens11 = ptens11 + rab(1)*fr(1) 474 ptens21 = ptens21 + rab(2)*fr(1) 475 ptens31 = ptens31 + rab(3)*fr(1) 476 ptens12 = ptens12 + rab(1)*fr(2) 477 ptens22 = ptens22 + rab(2)*fr(2) 478 ptens32 = ptens32 + rab(3)*fr(2) 479 ptens13 = ptens13 + rab(1)*fr(3) 480 ptens23 = ptens23 + rab(2)*fr(3) 481 ptens33 = ptens33 + rab(3)*fr(3) 482 END IF 483 ENDIF 484 END DO 485 END DO Kind_Group_Loop1 486 END DO 487 DEALLOCATE (eam_kinds_index) 488 489 ! Special way of handling the tersoff potential.. 490 IF (any_tersoff) THEN 491 NULLIFY (glob_loc_list, glob_cell_v, glob_loc_list_a) 492 CALL setup_tersoff_arrays(nonbonded, potparm, glob_loc_list, glob_cell_v, glob_loc_list_a, cell) 493 DO ilist = 1, nonbonded%nlists 494 neighbor_kind_pair => nonbonded%neighbor_kind_pairs(ilist) 495 npairs = neighbor_kind_pair%npairs 496 IF (npairs == 0) CYCLE 497 Kind_Group_Loop2: DO igrp = 1, neighbor_kind_pair%ngrp_kind 498 istart = neighbor_kind_pair%grp_kind_start(igrp) 499 iend = neighbor_kind_pair%grp_kind_end(igrp) 500 ikind = neighbor_kind_pair%ij_kind(1, igrp) 501 jkind = neighbor_kind_pair%ij_kind(2, igrp) 502 list => neighbor_kind_pair%list 503 cvi = neighbor_kind_pair%cell_vector 504 pot => potparm%pot(ikind, jkind)%pot 505 506 IF (pot%no_mb) CYCLE 507 rab2_max = pot%rcutsq 508 CALL matvec_3x3(cell_v, cell%hmat, cvi) 509 DO i = 1, SIZE(pot%type) 510 ! TERSOFF 511 IF (pot%type(i) == tersoff_type) THEN 512 npairs = iend - istart + 1 513 tersoff => pot%set(i)%tersoff 514 ALLOCATE (sort_list(2, npairs), work_list(npairs)) 515 sort_list = list(:, istart:iend) 516 ! Sort the list of neighbors, this increases the efficiency for single 517 ! potential contributions 518 CALL sort(sort_list(1, :), npairs, work_list) 519 DO ipair = 1, npairs 520 work_list(ipair) = sort_list(2, work_list(ipair)) 521 END DO 522 sort_list(2, :) = work_list 523 ! find number of unique elements of array index 1 524 nunique = 1 525 DO ipair = 1, npairs - 1 526 IF (sort_list(1, ipair + 1) /= sort_list(1, ipair)) nunique = nunique + 1 527 END DO 528 ipair = 1 529 junique = sort_list(1, ipair) 530 ifirst = 1 531 DO iunique = 1, nunique 532 atom_a = junique 533 IF (glob_loc_list_a(ifirst) > atom_a) CYCLE 534 DO mpair = ifirst, SIZE(glob_loc_list_a) 535 IF (glob_loc_list_a(mpair) == atom_a) EXIT 536 END DO 537 ifirst = mpair 538 DO mpair = ifirst, SIZE(glob_loc_list_a) 539 IF (glob_loc_list_a(mpair) /= atom_a) EXIT 540 END DO 541 ilast = mpair - 1 542 nloc_size = 0 543 IF (ifirst /= 0) nloc_size = ilast - ifirst + 1 544 DO WHILE (ipair <= npairs) 545 IF (sort_list(1, ipair) /= junique) EXIT 546 atom_b = sort_list(2, ipair) 547 ! Derivative terms 548 rtmp = r_last_update_pbc(atom_b)%r(:) - r_last_update_pbc(atom_a)%r(:) + cell_v 549 ipair = ipair + 1 550 IF (DOT_PRODUCT(rtmp, rtmp) <= tersoff%rcutsq) THEN 551 CALL tersoff_forces(tersoff, r_last_update_pbc, cell_v, & 552 nloc_size, glob_loc_list(:, ifirst:ilast), glob_cell_v(:, ifirst:ilast), & 553 atom_a, atom_b, f_nonbond, pv_nonbond, use_virial, tersoff%rcutsq) 554 END IF 555 END DO 556 ifirst = ilast + 1 557 IF (ipair <= npairs) junique = sort_list(1, ipair) 558 END DO 559 DEALLOCATE (sort_list, work_list) 560 END IF 561 END DO 562 END DO Kind_Group_Loop2 563 END DO 564 CALL destroy_tersoff_arrays(glob_loc_list, glob_cell_v, glob_loc_list_a) 565 END IF 566 ! Special way of handling the siepmann potential.. 567 IF (any_siepmann) THEN 568 NULLIFY (glob_loc_list, glob_cell_v, glob_loc_list_a) 569 CALL setup_siepmann_arrays(nonbonded, potparm, glob_loc_list, glob_cell_v, glob_loc_list_a, cell) 570 DO ilist = 1, nonbonded%nlists 571 neighbor_kind_pair => nonbonded%neighbor_kind_pairs(ilist) 572 npairs = neighbor_kind_pair%npairs 573 IF (npairs == 0) CYCLE 574 Kind_Group_Loop3: DO igrp = 1, neighbor_kind_pair%ngrp_kind 575 istart = neighbor_kind_pair%grp_kind_start(igrp) 576 iend = neighbor_kind_pair%grp_kind_end(igrp) 577 ikind = neighbor_kind_pair%ij_kind(1, igrp) 578 jkind = neighbor_kind_pair%ij_kind(2, igrp) 579 list => neighbor_kind_pair%list 580 cvi = neighbor_kind_pair%cell_vector 581 pot => potparm%pot(ikind, jkind)%pot 582 583 IF (pot%no_mb) CYCLE 584 rab2_max = pot%rcutsq 585 CALL matvec_3x3(cell_v, cell%hmat, cvi) 586 DO i = 1, SIZE(pot%type) 587 ! SIEPMANN 588 IF (pot%type(i) == siepmann_type) THEN 589 npairs = iend - istart + 1 590 siepmann => pot%set(i)%siepmann 591 ALLOCATE (sort_list(2, npairs), work_list(npairs)) 592 sort_list = list(:, istart:iend) 593 ! Sort the list of neighbors, this increases the efficiency for single 594 ! potential contributions 595 CALL sort(sort_list(1, :), npairs, work_list) 596 DO ipair = 1, npairs 597 work_list(ipair) = sort_list(2, work_list(ipair)) 598 END DO 599 sort_list(2, :) = work_list 600 ! find number of unique elements of array index 1 601 nunique = 1 602 DO ipair = 1, npairs - 1 603 IF (sort_list(1, ipair + 1) /= sort_list(1, ipair)) nunique = nunique + 1 604 END DO 605 ipair = 1 606 junique = sort_list(1, ipair) 607 ifirst = 1 608 DO iunique = 1, nunique 609 atom_a = junique 610 IF (glob_loc_list_a(ifirst) > atom_a) CYCLE 611 DO mpair = ifirst, SIZE(glob_loc_list_a) 612 IF (glob_loc_list_a(mpair) == atom_a) EXIT 613 END DO 614 ifirst = mpair 615 DO mpair = ifirst, SIZE(glob_loc_list_a) 616 IF (glob_loc_list_a(mpair) /= atom_a) EXIT 617 END DO 618 ilast = mpair - 1 619 nloc_size = 0 620 IF (ifirst /= 0) nloc_size = ilast - ifirst + 1 621 DO WHILE (ipair <= npairs) 622 IF (sort_list(1, ipair) /= junique) EXIT 623 atom_b = sort_list(2, ipair) 624 ! Derivative terms 625 rtmp = r_last_update_pbc(atom_b)%r(:) - r_last_update_pbc(atom_a)%r(:) + cell_v 626 ipair = ipair + 1 627 IF (DOT_PRODUCT(rtmp, rtmp) <= siepmann%rcutsq) THEN 628 CALL siepmann_forces_v2(siepmann, r_last_update_pbc, cell_v, cell, & 629 atom_a, atom_b, f_nonbond, use_virial, siepmann%rcutsq, & 630 particle_set) 631 CALL siepmann_forces_v3(siepmann, r_last_update_pbc, cell_v, & 632 nloc_size, glob_loc_list(:, ifirst:ilast), & 633 atom_a, atom_b, f_nonbond, use_virial, siepmann%rcutsq, & 634 cell, particle_set) 635 END IF 636 END DO 637 ifirst = ilast + 1 638 IF (ipair <= npairs) junique = sort_list(1, ipair) 639 END DO 640 DEALLOCATE (sort_list, work_list) 641 END IF 642 END DO 643 END DO Kind_Group_Loop3 644 END DO 645 CALL destroy_siepmann_arrays(glob_loc_list, glob_cell_v, glob_loc_list_a) 646 END IF 647 IF (use_virial) THEN 648 pv_nonbond(1, 1) = pv_nonbond(1, 1) + ptens11 649 pv_nonbond(1, 2) = pv_nonbond(1, 2) + ptens12 650 pv_nonbond(1, 3) = pv_nonbond(1, 3) + ptens13 651 pv_nonbond(2, 1) = pv_nonbond(2, 1) + ptens21 652 pv_nonbond(2, 2) = pv_nonbond(2, 2) + ptens22 653 pv_nonbond(2, 3) = pv_nonbond(2, 3) + ptens23 654 pv_nonbond(3, 1) = pv_nonbond(3, 1) + ptens31 655 pv_nonbond(3, 2) = pv_nonbond(3, 2) + ptens32 656 pv_nonbond(3, 3) = pv_nonbond(3, 3) + ptens33 657 END IF 658 CALL timestop(handle) 659 END SUBROUTINE force_nonbond_manybody 660 661END MODULE manybody_potential 662 663