|
CP2K 2.4 (Revision 12889)
|
00001 !-----------------------------------------------------------------------------! 00002 ! CP2K: A general program to perform molecular dynamics simulations ! 00003 ! Copyright (C) 2000 - 2013 CP2K developers group ! 00004 !-----------------------------------------------------------------------------! 00005 00006 ! ***************************************************************************** 00014 MODULE manybody_potential 00015 00016 USE atomic_kind_types, ONLY: atomic_kind_type 00017 USE atprop_types, ONLY: atprop_type 00018 USE cell_types, ONLY: cell_type 00019 USE cp_para_types, ONLY: cp_para_env_type 00020 USE distribution_1d_types, ONLY: distribution_1d_type 00021 USE f77_blas 00022 USE fist_neighbor_list_types, ONLY: fist_neighbor_type,& 00023 neighbor_kind_pairs_type 00024 USE fist_nonbond_env_types, ONLY: eam_type,& 00025 fist_nonbond_env_get,& 00026 fist_nonbond_env_type,& 00027 pos_type 00028 USE kinds, ONLY: dp 00029 USE manybody_eam, ONLY: get_force_eam 00030 USE manybody_siepmann, ONLY: destroy_siepmann_arrays,& 00031 setup_siepmann_arrays,& 00032 siepmann_energy,& 00033 siepmann_forces_v2,& 00034 siepmann_forces_v3 00035 USE manybody_tersoff, ONLY: destroy_tersoff_arrays,& 00036 setup_tersoff_arrays,& 00037 tersoff_energy,& 00038 tersoff_forces 00039 USE mathlib, ONLY: matvec_3x3 00040 USE message_passing, ONLY: mp_sum 00041 USE pair_potential_types, ONLY: ea_type,& 00042 eam_pot_type,& 00043 pair_potential_pp_type,& 00044 pair_potential_single_type,& 00045 siepmann_pot_type,& 00046 siepmann_type,& 00047 tersoff_pot_type,& 00048 tersoff_type 00049 USE particle_types, ONLY: particle_type 00050 USE timings, ONLY: timeset,& 00051 timestop 00052 USE util, ONLY: sort 00053 #include "cp_common_uses.h" 00054 00055 IMPLICIT NONE 00056 00057 PRIVATE 00058 PUBLIC :: energy_manybody 00059 PUBLIC :: force_nonbond_manybody 00060 CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'manybody_potential' 00061 00062 CONTAINS 00063 00064 ! ***************************************************************************** 00072 SUBROUTINE energy_manybody ( fist_nonbond_env, atomic_kind_set, local_particles,& 00073 particle_set, cell, pot_manybody, atprop_env, para_env, error ) 00074 00075 TYPE(fist_nonbond_env_type), POINTER :: fist_nonbond_env 00076 TYPE(atomic_kind_type), POINTER :: atomic_kind_set( : ) 00077 TYPE(distribution_1d_type), POINTER :: local_particles 00078 TYPE(particle_type), POINTER :: particle_set( : ) 00079 TYPE(cell_type), POINTER :: cell 00080 REAL(dp), INTENT(INOUT) :: pot_manybody 00081 TYPE(atprop_type), POINTER :: atprop_env 00082 TYPE(cp_para_env_type), OPTIONAL, 00083 POINTER :: para_env 00084 TYPE(cp_error_type), INTENT(inout) :: error 00085 00086 CHARACTER(LEN=*), PARAMETER :: routineN = 'energy_manybody', 00087 routineP = moduleN//':'//routineN 00088 00089 INTEGER :: atom_a, atom_b, handle, i, iend, ifirst, igrp, ikind, ilast, 00090 ilist, indexa, ipair, iparticle, iparticle_local, istart, iunique, 00091 jkind, junique, mpair, nkinds, nloc_size, npairs, nparticle, 00092 nparticle_local, nunique, stat 00093 INTEGER, DIMENSION(:), POINTER :: glob_loc_list_a, work_list 00094 INTEGER, DIMENSION(:, :), POINTER :: glob_loc_list, list, sort_list 00095 LOGICAL :: any_siepmann, any_tersoff, 00096 failure 00097 REAL(KIND=dp) :: drij, embed, pot_loc, qr, 00098 rab2_max, rij(3) 00099 REAL(KIND=dp), DIMENSION(3) :: cell_v, cvi 00100 REAL(KIND=dp), DIMENSION(:, :), POINTER :: glob_cell_v 00101 REAL(KIND=dp), POINTER :: fembed( : ) 00102 TYPE(eam_pot_type), POINTER :: eam 00103 TYPE(eam_type), DIMENSION(:), POINTER :: eam_data 00104 TYPE(fist_neighbor_type), POINTER :: nonbonded 00105 TYPE(neighbor_kind_pairs_type), POINTER :: neighbor_kind_pair 00106 TYPE(pair_potential_pp_type), POINTER :: potparm 00107 TYPE(pair_potential_single_type), 00108 POINTER :: pot 00109 TYPE(pos_type), DIMENSION(:), POINTER :: r_last_update_pbc 00110 TYPE(siepmann_pot_type), POINTER :: siepmann 00111 TYPE(tersoff_pot_type), POINTER :: tersoff 00112 00113 NULLIFY(eam, siepmann, tersoff) 00114 any_tersoff = .FALSE. 00115 any_siepmann = .FALSE. 00116 failure = .FALSE. 00117 CALL timeset ( routineN, handle ) 00118 CALL fist_nonbond_env_get ( fist_nonbond_env, r_last_update_pbc=r_last_update_pbc,& 00119 potparm = potparm , eam_data=eam_data, error=error) 00120 ! EAM requires a single loop 00121 DO ikind = 1, SIZE ( atomic_kind_set ) 00122 pot => potparm %pot ( ikind, ikind ) % pot 00123 DO i = 1, SIZE(pot%type) 00124 IF(pot%type(i) /= ea_type) CYCLE 00125 eam => pot%set(i)%eam 00126 nparticle = SIZE ( particle_set ) 00127 ALLOCATE(fembed(nparticle), stat=stat ) 00128 CPPostcondition(stat==0,cp_failure_level,routineP,error,failure) 00129 fembed(:) = 0._dp 00130 CPPostcondition(ASSOCIATED(eam_data),cp_failure_level,routineP,error,failure) 00131 ! computation of embedding function and energy 00132 nparticle_local = local_particles%n_el(ikind) 00133 DO iparticle_local=1,nparticle_local 00134 iparticle = local_particles%list(ikind)%array(iparticle_local) 00135 indexa = INT ( eam_data(iparticle)%rho / eam%drhoar) + 1 00136 IF (indexa > eam%npoints-1) indexa = eam%npoints-1 00137 qr = eam_data(iparticle)%rho - eam%rhoval(indexa) 00138 00139 embed = eam%frho(indexa) + qr*eam%frhop(indexa) 00140 fembed (iparticle) = eam%frhop(indexa) + qr*(eam%frhop(indexa+1)-eam%frhop(indexa))/eam%drhoar 00141 00142 pot_manybody = pot_manybody + embed 00143 END DO 00144 ! communicate data 00145 CALL mp_sum(fembed, para_env%group) 00146 DO iparticle=1,nparticle 00147 IF (particle_set(iparticle)%atomic_kind%kind_number==ikind) THEN 00148 eam_data(iparticle)%f_embed = fembed(iparticle) 00149 END IF 00150 END DO 00151 00152 DEALLOCATE ( fembed, stat=stat ) 00153 CPPostcondition(stat==0,cp_failure_level,routineP,error,failure) 00154 END DO 00155 END DO 00156 ! Other manybody potential 00157 DO ikind = 1, SIZE ( atomic_kind_set ) 00158 DO jkind = ikind, SIZE ( atomic_kind_set ) 00159 pot => potparm %pot ( ikind, jkind ) % pot 00160 any_tersoff = any_tersoff .OR. ANY(pot%type==tersoff_type) 00161 any_siepmann = any_siepmann .OR. ANY(pot%type==siepmann_type) 00162 END DO 00163 END DO 00164 CALL fist_nonbond_env_get(fist_nonbond_env,nonbonded=nonbonded,natom_types=nkinds,error=error) 00165 ! TERSOFF 00166 IF (any_tersoff) THEN 00167 NULLIFY(glob_loc_list, glob_cell_v, glob_loc_list_a) 00168 CALL setup_tersoff_arrays(nonbonded, potparm, glob_loc_list, glob_cell_v, glob_loc_list_a, cell, error) 00169 DO ilist=1,nonbonded%nlists 00170 neighbor_kind_pair => nonbonded%neighbor_kind_pairs(ilist) 00171 npairs=neighbor_kind_pair%npairs 00172 IF (npairs ==0) CYCLE 00173 Kind_Group_Loop: DO igrp = 1, neighbor_kind_pair%ngrp_kind 00174 istart = neighbor_kind_pair%grp_kind_start(igrp) 00175 iend = neighbor_kind_pair%grp_kind_end(igrp) 00176 ikind = neighbor_kind_pair%ij_kind(1,igrp) 00177 jkind = neighbor_kind_pair%ij_kind(2,igrp) 00178 list => neighbor_kind_pair%list 00179 cvi = neighbor_kind_pair%cell_vector 00180 pot => potparm %pot ( ikind, jkind ) % pot 00181 DO i = 1, SIZE(pot%type) 00182 IF(pot%type(i) /= tersoff_type) CYCLE 00183 rab2_max = pot%set(i)%tersoff%rcutsq 00184 CALL matvec_3x3(cell_v, cell%hmat,cvi) 00185 pot => potparm%pot( ikind, jkind )%pot 00186 tersoff => pot%set(i)%tersoff 00187 npairs = iend-istart+1 00188 IF (npairs /=0) THEN 00189 ALLOCATE(sort_list(2,npairs),work_list (npairs),stat=stat) 00190 CPPostcondition(stat==0,cp_failure_level,routineP,error,failure) 00191 sort_list = list(:,istart:iend) 00192 ! Sort the list of neighbors, this increases the efficiency for single 00193 ! potential contributions 00194 CALL sort(sort_list(1,:),npairs,work_list) 00195 DO ipair = 1, npairs 00196 work_list(ipair)=sort_list(2,work_list(ipair)) 00197 END DO 00198 sort_list(2,:) = work_list 00199 ! find number of unique elements of array index 1 00200 nunique = 1 00201 DO ipair = 1, npairs-1 00202 IF (sort_list(1,ipair+1)/=sort_list(1,ipair)) nunique = nunique + 1 00203 END DO 00204 ipair = 1 00205 junique = sort_list(1,ipair) 00206 ifirst = 1 00207 DO iunique = 1, nunique 00208 atom_a = junique 00209 IF (glob_loc_list_a(ifirst)>atom_a) CYCLE 00210 DO mpair = ifirst, SIZE(glob_loc_list_a) 00211 IF (glob_loc_list_a(mpair)==atom_a) EXIT 00212 END DO 00213 ifirst = mpair 00214 DO mpair = ifirst, SIZE(glob_loc_list_a) 00215 IF (glob_loc_list_a(mpair)/=atom_a) EXIT 00216 END DO 00217 ilast = mpair-1 00218 nloc_size = 0 00219 IF (ifirst/=0) nloc_size = ilast-ifirst+1 00220 DO WHILE (ipair<=npairs) 00221 IF (sort_list(1,ipair)/=junique) EXIT 00222 atom_b = sort_list(2,ipair) 00223 ! Energy terms 00224 pot_loc = 0.0_dp 00225 rij(:) = r_last_update_pbc(atom_b)%r(:)-r_last_update_pbc(atom_a)%r(:)+cell_v 00226 drij = DOT_PRODUCT(rij,rij) 00227 ipair = ipair + 1 00228 IF(drij>rab2_max) CYCLE 00229 drij = SQRT(drij) 00230 CALL tersoff_energy(pot_loc, tersoff, r_last_update_pbc, atom_a, atom_b, nloc_size,& 00231 glob_loc_list(:,ifirst:ilast),glob_cell_v(:,ifirst:ilast), cell_v, drij) 00232 pot_manybody = pot_manybody + 0.5_dp*pot_loc 00233 END DO 00234 ifirst = ilast + 1 00235 IF (ipair<=npairs) junique = sort_list(1,ipair) 00236 END DO 00237 DEALLOCATE(sort_list,work_list,stat=stat) 00238 CPPostcondition(stat==0,cp_failure_level,routineP,error,failure) 00239 END IF 00240 END DO 00241 END DO Kind_Group_Loop 00242 END DO 00243 CALL destroy_tersoff_arrays(glob_loc_list, glob_cell_v, glob_loc_list_a, error) 00244 END IF 00245 00246 !SIEPMANN POTENTIAL 00247 IF (any_siepmann) THEN 00248 NULLIFY(glob_loc_list, glob_cell_v, glob_loc_list_a) 00249 CALL setup_siepmann_arrays(nonbonded, potparm, glob_loc_list, glob_cell_v, glob_loc_list_a, cell, error) 00250 DO ilist=1,nonbonded%nlists 00251 neighbor_kind_pair => nonbonded%neighbor_kind_pairs(ilist) 00252 npairs=neighbor_kind_pair%npairs 00253 IF (npairs ==0) CYCLE 00254 Kind_Group_Loop_2: DO igrp = 1, neighbor_kind_pair%ngrp_kind 00255 istart = neighbor_kind_pair%grp_kind_start(igrp) 00256 iend = neighbor_kind_pair%grp_kind_end(igrp) 00257 ikind = neighbor_kind_pair%ij_kind(1,igrp) 00258 jkind = neighbor_kind_pair%ij_kind(2,igrp) 00259 list => neighbor_kind_pair%list 00260 cvi = neighbor_kind_pair%cell_vector 00261 pot => potparm %pot ( ikind, jkind ) % pot 00262 DO i = 1, SIZE(pot%type) 00263 IF(pot%type(i) /= siepmann_type) CYCLE 00264 rab2_max = pot%set(i)%siepmann%rcutsq 00265 CALL matvec_3x3(cell_v, cell%hmat,cvi) 00266 pot => potparm%pot( ikind, jkind )%pot 00267 siepmann => pot%set(i)%siepmann 00268 npairs = iend-istart+1 00269 IF (npairs /=0) THEN 00270 ALLOCATE(sort_list(2,npairs),work_list (npairs),stat=stat) 00271 CPPostcondition(stat==0,cp_failure_level,routineP,error,failure) 00272 sort_list = list(:,istart:iend) 00273 ! Sort the list of neighbors, this increases the efficiency for single 00274 ! potential contributions 00275 CALL sort(sort_list(1,:),npairs,work_list) 00276 DO ipair = 1, npairs 00277 work_list(ipair)=sort_list(2,work_list(ipair)) 00278 END DO 00279 sort_list(2,:) = work_list 00280 ! find number of unique elements of array index 1 00281 nunique = 1 00282 DO ipair = 1, npairs-1 00283 IF (sort_list(1,ipair+1)/=sort_list(1,ipair)) nunique = nunique + 1 00284 END DO 00285 ipair = 1 00286 junique = sort_list(1,ipair) 00287 ifirst = 1 00288 DO iunique = 1, nunique 00289 atom_a = junique 00290 IF (glob_loc_list_a(ifirst)>atom_a) CYCLE 00291 DO mpair = ifirst, SIZE(glob_loc_list_a) 00292 IF (glob_loc_list_a(mpair)==atom_a) EXIT 00293 END DO 00294 ifirst = mpair 00295 DO mpair = ifirst, SIZE(glob_loc_list_a) 00296 IF (glob_loc_list_a(mpair)/=atom_a) EXIT 00297 END DO 00298 ilast = mpair-1 00299 nloc_size = 0 00300 IF (ifirst/=0) nloc_size = ilast-ifirst+1 00301 DO WHILE (ipair<=npairs) 00302 IF (sort_list(1,ipair)/=junique) EXIT 00303 atom_b = sort_list(2,ipair) 00304 ! Energy terms 00305 pot_loc = 0.0_dp 00306 rij(:) = r_last_update_pbc(atom_b)%r(:)-r_last_update_pbc(atom_a)%r(:)+cell_v 00307 drij = DOT_PRODUCT(rij,rij) 00308 ipair = ipair + 1 00309 IF(drij>rab2_max) CYCLE 00310 drij = SQRT(drij) 00311 CALL siepmann_energy(pot_loc, siepmann, r_last_update_pbc, atom_a, atom_b, nloc_size,& 00312 glob_loc_list(:,ifirst:ilast),glob_cell_v(:,ifirst:ilast), cell_v, cell, drij, particle_set) 00313 pot_manybody = pot_manybody + pot_loc 00314 END DO 00315 ifirst = ilast + 1 00316 IF (ipair<=npairs) junique = sort_list(1,ipair) 00317 END DO 00318 DEALLOCATE(sort_list,work_list,stat=stat) 00319 CPPostcondition(stat==0,cp_failure_level,routineP,error,failure) 00320 END IF 00321 END DO 00322 END DO Kind_Group_Loop_2 00323 END DO 00324 CALL destroy_siepmann_arrays(glob_loc_list, glob_cell_v, glob_loc_list_a, error) 00325 END IF 00326 CALL timestop ( handle ) 00327 END SUBROUTINE energy_manybody 00328 00329 ! ***************************************************************************** 00336 SUBROUTINE force_nonbond_manybody ( fist_nonbond_env, particle_set, cell, & 00337 f_nonbond, pv_nonbond, atprop_env, use_virial, error ) 00338 00339 TYPE(fist_nonbond_env_type), POINTER :: fist_nonbond_env 00340 TYPE(particle_type), DIMENSION(:), 00341 POINTER :: particle_set 00342 TYPE(cell_type), POINTER :: cell 00343 REAL(KIND=dp), DIMENSION(:, :), 00344 INTENT(INOUT) :: f_nonbond, pv_nonbond 00345 TYPE(atprop_type), POINTER :: atprop_env 00346 LOGICAL, INTENT(IN) :: use_virial 00347 TYPE(cp_error_type), INTENT(inout) :: error 00348 00349 CHARACTER(LEN=*), PARAMETER :: routineN = 'force_nonbond_manybody', 00350 routineP = moduleN//':'//routineN 00351 00352 INTEGER :: atom_a, atom_b, handle, i, i_a, i_b, iend, ifirst, igrp, 00353 ikind, ilast, ilist, ipair, istart, iunique, jkind, junique, kind_a, 00354 kind_b, mpair, nkinds, nloc_size, npairs, nunique, stat 00355 INTEGER, ALLOCATABLE, DIMENSION(:, :) :: eam_kinds_index 00356 INTEGER, DIMENSION(:), POINTER :: glob_loc_list_a, work_list 00357 INTEGER, DIMENSION(:, :), POINTER :: glob_loc_list, list, sort_list 00358 LOGICAL :: any_siepmann, any_tersoff, 00359 failure 00360 REAL(KIND=dp) :: f_eam, fac, fr(3), ptens11, ptens12, ptens13, ptens21, 00361 ptens22, ptens23, ptens31, ptens32, ptens33, rab(3), rab2, rab2_max, 00362 rtmp(3) 00363 REAL(KIND=dp), DIMENSION(3) :: cell_v, cvi 00364 REAL(KIND=dp), DIMENSION(:, :), POINTER :: glob_cell_v 00365 TYPE(eam_pot_type), POINTER :: eam_a, eam_b 00366 TYPE(eam_type), DIMENSION(:), POINTER :: eam_data 00367 TYPE(fist_neighbor_type), POINTER :: nonbonded 00368 TYPE(neighbor_kind_pairs_type), POINTER :: neighbor_kind_pair 00369 TYPE(pair_potential_pp_type), POINTER :: potparm 00370 TYPE(pair_potential_single_type), 00371 POINTER :: pot 00372 TYPE(pos_type), DIMENSION(:), POINTER :: r_last_update_pbc 00373 TYPE(siepmann_pot_type), POINTER :: siepmann 00374 TYPE(tersoff_pot_type), POINTER :: tersoff 00375 00376 failure = .FALSE. 00377 any_tersoff = .FALSE. 00378 any_siepmann = .FALSE. 00379 CALL timeset ( routineN, handle ) 00380 NULLIFY(eam_a, eam_b, tersoff, siepmann) 00381 00382 CALL fist_nonbond_env_get(fist_nonbond_env,nonbonded=nonbonded,potparm=potparm,& 00383 natom_types=nkinds, eam_data=eam_data, r_last_update_pbc=r_last_update_pbc,error=error) 00384 00385 ! Initializing the potential energy, pressure tensor and force 00386 IF (use_virial) THEN 00387 ptens11 = 0.0_dp ; ptens12 = 0.0_dp ; ptens13 = 0.0_dp 00388 ptens21 = 0.0_dp ; ptens22 = 0.0_dp ; ptens23 = 0.0_dp 00389 ptens31 = 0.0_dp ; ptens32 = 0.0_dp ; ptens33 = 0.0_dp 00390 END IF 00391 00392 nkinds = SIZE(potparm%pot,1) 00393 ALLOCATE(eam_kinds_index(nkinds,nkinds),stat=stat) 00394 CPPostcondition(stat==0,cp_failure_level,routineP,error,failure) 00395 eam_kinds_index = -1 00396 DO ikind = 1, nkinds 00397 DO jkind = ikind, nkinds 00398 DO i = 1, SIZE(potparm%pot(ikind,jkind)%pot%type) 00399 IF (potparm%pot(ikind,jkind)%pot%type(i)==ea_type) THEN 00400 ! At the moment we allow only 1 EAM per each kinds pair.. 00401 CPPostcondition(eam_kinds_index(ikind,jkind)==-1,cp_failure_level,routineP,error,failure) 00402 CPPostcondition(eam_kinds_index(jkind,ikind)==-1,cp_failure_level,routineP,error,failure) 00403 eam_kinds_index(ikind,jkind) = i 00404 eam_kinds_index(jkind,ikind) = i 00405 END IF 00406 END DO 00407 END DO 00408 END DO 00409 00410 ! starting the force loop 00411 DO ilist=1,nonbonded%nlists 00412 neighbor_kind_pair => nonbonded%neighbor_kind_pairs(ilist) 00413 npairs=neighbor_kind_pair%npairs 00414 IF (npairs ==0) CYCLE 00415 Kind_Group_Loop1: DO igrp = 1, neighbor_kind_pair%ngrp_kind 00416 istart = neighbor_kind_pair%grp_kind_start(igrp) 00417 iend = neighbor_kind_pair%grp_kind_end(igrp) 00418 ikind = neighbor_kind_pair%ij_kind(1,igrp) 00419 jkind = neighbor_kind_pair%ij_kind(2,igrp) 00420 list => neighbor_kind_pair%list 00421 cvi = neighbor_kind_pair%cell_vector 00422 pot => potparm %pot( ikind, jkind )%pot 00423 IF (pot%no_mb) CYCLE 00424 rab2_max = pot%rcutsq 00425 CALL matvec_3x3(cell_v, cell%hmat,cvi) 00426 any_tersoff = any_tersoff.OR.ANY(pot%type==tersoff_type) 00427 any_siepmann = any_siepmann.OR.ANY(pot%type==siepmann_type) 00428 i = eam_kinds_index(ikind,jkind) 00429 IF (i==-1) CYCLE 00430 ! EAM 00431 CPPostcondition(ASSOCIATED(eam_data),cp_failure_level,routineP,error,failure) 00432 DO ipair = istart, iend 00433 atom_a = list(1,ipair) 00434 atom_b = list(2,ipair) 00435 fac = 1.0_dp 00436 IF (atom_a==atom_b) fac = 0.5_dp 00437 kind_a = particle_set(atom_a)%atomic_kind%kind_number 00438 kind_b = particle_set(atom_b)%atomic_kind%kind_number 00439 i_a = eam_kinds_index(kind_a,kind_a) 00440 i_b = eam_kinds_index(kind_b,kind_b) 00441 eam_a => potparm%pot(kind_a,kind_a)%pot%set(i_a)%eam 00442 eam_b => potparm%pot(kind_b,kind_b)%pot%set(i_b)%eam 00443 00444 !set this outside the potential type in case need multiple potentials 00445 !Do everything necessary for EAM here 00446 rab = r_last_update_pbc(atom_b)%r-r_last_update_pbc(atom_a)%r 00447 rab = rab+cell_v 00448 rab2 = rab(1)*rab(1) + rab(2)*rab(2) + rab(3)*rab(3) 00449 IF (rab2 <= rab2_max) THEN 00450 CALL get_force_eam ( rab2, eam_a, eam_b, eam_data, atom_a, atom_b, f_eam ) 00451 f_eam = f_eam * fac 00452 00453 fr(1) = - f_eam * rab(1) 00454 fr(2) = - f_eam * rab(2) 00455 fr(3) = - f_eam * rab(3) 00456 f_nonbond(1,atom_a) = f_nonbond(1,atom_a) - fr(1) 00457 f_nonbond(2,atom_a) = f_nonbond(2,atom_a) - fr(2) 00458 f_nonbond(3,atom_a) = f_nonbond(3,atom_a) - fr(3) 00459 00460 f_nonbond(1,atom_b) = f_nonbond(1,atom_b) + fr(1) 00461 f_nonbond(2,atom_b) = f_nonbond(2,atom_b) + fr(2) 00462 f_nonbond(3,atom_b) = f_nonbond(3,atom_b) + fr(3) 00463 IF (use_virial) THEN 00464 ptens11 = ptens11 + rab(1)*fr(1) 00465 ptens21 = ptens21 + rab(2)*fr(1) 00466 ptens31 = ptens31 + rab(3)*fr(1) 00467 ptens12 = ptens12 + rab(1)*fr(2) 00468 ptens22 = ptens22 + rab(2)*fr(2) 00469 ptens32 = ptens32 + rab(3)*fr(2) 00470 ptens13 = ptens13 + rab(1)*fr(3) 00471 ptens23 = ptens23 + rab(2)*fr(3) 00472 ptens33 = ptens33 + rab(3)*fr(3) 00473 END IF 00474 ENDIF 00475 END DO 00476 END DO Kind_Group_Loop1 00477 END DO 00478 DEALLOCATE (eam_kinds_index, stat=stat) 00479 CPPostcondition(stat==0,cp_failure_level,routineP,error,failure) 00480 00481 ! Special way of handling the tersoff potential.. 00482 IF (any_tersoff) THEN 00483 NULLIFY(glob_loc_list, glob_cell_v, glob_loc_list_a) 00484 CALL setup_tersoff_arrays(nonbonded, potparm, glob_loc_list, glob_cell_v, glob_loc_list_a, cell, error) 00485 DO ilist=1,nonbonded%nlists 00486 neighbor_kind_pair => nonbonded%neighbor_kind_pairs(ilist) 00487 npairs=neighbor_kind_pair%npairs 00488 IF (npairs ==0) CYCLE 00489 Kind_Group_Loop2: DO igrp = 1, neighbor_kind_pair%ngrp_kind 00490 istart = neighbor_kind_pair%grp_kind_start(igrp) 00491 iend = neighbor_kind_pair%grp_kind_end(igrp) 00492 ikind = neighbor_kind_pair%ij_kind(1,igrp) 00493 jkind = neighbor_kind_pair%ij_kind(2,igrp) 00494 list => neighbor_kind_pair%list 00495 cvi = neighbor_kind_pair%cell_vector 00496 pot => potparm %pot( ikind, jkind )%pot 00497 00498 IF (pot%no_mb) CYCLE 00499 rab2_max = pot%rcutsq 00500 CALL matvec_3x3(cell_v, cell%hmat,cvi) 00501 DO i = 1, SIZE(pot%type) 00502 ! TERSOFF 00503 IF (pot%type(i)==tersoff_type) THEN 00504 npairs = iend-istart+1 00505 tersoff => pot%set(i)%tersoff 00506 ALLOCATE(sort_list(2,npairs),work_list (npairs),stat=stat) 00507 CPPostcondition(stat==0,cp_failure_level,routineP,error,failure) 00508 sort_list = list(:,istart:iend) 00509 ! Sort the list of neighbors, this increases the efficiency for single 00510 ! potential contributions 00511 CALL sort(sort_list(1,:),npairs,work_list) 00512 DO ipair = 1, npairs 00513 work_list(ipair)=sort_list(2,work_list(ipair)) 00514 END DO 00515 sort_list(2,:) = work_list 00516 ! find number of unique elements of array index 1 00517 nunique = 1 00518 DO ipair = 1, npairs-1 00519 IF (sort_list(1,ipair+1)/=sort_list(1,ipair)) nunique = nunique + 1 00520 END DO 00521 ipair = 1 00522 junique = sort_list(1,ipair) 00523 ifirst = 1 00524 DO iunique = 1, nunique 00525 atom_a = junique 00526 IF (glob_loc_list_a(ifirst)>atom_a) CYCLE 00527 DO mpair = ifirst, SIZE(glob_loc_list_a) 00528 IF (glob_loc_list_a(mpair)==atom_a) EXIT 00529 END DO 00530 ifirst = mpair 00531 DO mpair = ifirst, SIZE(glob_loc_list_a) 00532 IF (glob_loc_list_a(mpair)/=atom_a) EXIT 00533 END DO 00534 ilast = mpair-1 00535 nloc_size = 0 00536 IF (ifirst/=0) nloc_size = ilast-ifirst+1 00537 DO WHILE (ipair<=npairs) 00538 IF (sort_list(1,ipair)/=junique) EXIT 00539 atom_b = sort_list(2,ipair) 00540 ! Derivative terms 00541 rtmp = r_last_update_pbc(atom_b)%r(:)-r_last_update_pbc(atom_a)%r(:) + cell_v 00542 ipair = ipair + 1 00543 IF (DOT_PRODUCT(rtmp,rtmp)<=tersoff%rcutsq) THEN 00544 CALL tersoff_forces ( tersoff, r_last_update_pbc, cell_v,& 00545 nloc_size, glob_loc_list(:,ifirst:ilast), glob_cell_v(:,ifirst:ilast),& 00546 atom_a, atom_b, f_nonbond, pv_nonbond, use_virial, tersoff%rcutsq) 00547 END IF 00548 END DO 00549 ifirst = ilast + 1 00550 IF (ipair<=npairs) junique = sort_list(1,ipair) 00551 END DO 00552 DEALLOCATE(sort_list,work_list,stat=stat) 00553 CPPostcondition(stat==0,cp_failure_level,routineP,error,failure) 00554 END IF 00555 END DO 00556 END DO Kind_Group_Loop2 00557 END DO 00558 CALL destroy_tersoff_arrays(glob_loc_list, glob_cell_v, glob_loc_list_a, error) 00559 END IF 00560 ! Special way of handling the siepmann potential.. 00561 IF (any_siepmann) THEN 00562 NULLIFY(glob_loc_list, glob_cell_v, glob_loc_list_a) 00563 CALL setup_siepmann_arrays(nonbonded, potparm, glob_loc_list, glob_cell_v, glob_loc_list_a, cell, error) 00564 DO ilist=1,nonbonded%nlists 00565 neighbor_kind_pair => nonbonded%neighbor_kind_pairs(ilist) 00566 npairs=neighbor_kind_pair%npairs 00567 IF (npairs ==0) CYCLE 00568 Kind_Group_Loop3: DO igrp = 1, neighbor_kind_pair%ngrp_kind 00569 istart = neighbor_kind_pair%grp_kind_start(igrp) 00570 iend = neighbor_kind_pair%grp_kind_end(igrp) 00571 ikind = neighbor_kind_pair%ij_kind(1,igrp) 00572 jkind = neighbor_kind_pair%ij_kind(2,igrp) 00573 list => neighbor_kind_pair%list 00574 cvi = neighbor_kind_pair%cell_vector 00575 pot => potparm %pot( ikind, jkind )%pot 00576 00577 IF (pot%no_mb) CYCLE 00578 rab2_max = pot%rcutsq 00579 CALL matvec_3x3(cell_v, cell%hmat,cvi) 00580 DO i = 1, SIZE(pot%type) 00581 ! SIEPMANN 00582 IF (pot%type(i)==siepmann_type) THEN 00583 npairs = iend-istart+1 00584 siepmann => pot%set(i)%siepmann 00585 ALLOCATE(sort_list(2,npairs),work_list (npairs),stat=stat) 00586 CPPostcondition(stat==0,cp_failure_level,routineP,error,failure) 00587 sort_list = list(:,istart:iend) 00588 ! Sort the list of neighbors, this increases the efficiency for single 00589 ! potential contributions 00590 CALL sort(sort_list(1,:),npairs,work_list) 00591 DO ipair = 1, npairs 00592 work_list(ipair)=sort_list(2,work_list(ipair)) 00593 END DO 00594 sort_list(2,:) = work_list 00595 ! find number of unique elements of array index 1 00596 nunique = 1 00597 DO ipair = 1, npairs-1 00598 IF (sort_list(1,ipair+1)/=sort_list(1,ipair)) nunique = nunique + 1 00599 END DO 00600 ipair = 1 00601 junique = sort_list(1,ipair) 00602 ifirst = 1 00603 DO iunique = 1, nunique 00604 atom_a = junique 00605 IF (glob_loc_list_a(ifirst)>atom_a) CYCLE 00606 DO mpair = ifirst, SIZE(glob_loc_list_a) 00607 IF (glob_loc_list_a(mpair)==atom_a) EXIT 00608 END DO 00609 ifirst = mpair 00610 DO mpair = ifirst, SIZE(glob_loc_list_a) 00611 IF (glob_loc_list_a(mpair)/=atom_a) EXIT 00612 END DO 00613 ilast = mpair-1 00614 nloc_size = 0 00615 IF (ifirst/=0) nloc_size = ilast-ifirst+1 00616 DO WHILE (ipair<=npairs) 00617 IF (sort_list(1,ipair)/=junique) EXIT 00618 atom_b = sort_list(2,ipair) 00619 ! Derivative terms 00620 rtmp = r_last_update_pbc(atom_b)%r(:)-r_last_update_pbc(atom_a)%r(:) + cell_v 00621 ipair = ipair + 1 00622 IF (DOT_PRODUCT(rtmp,rtmp)<=siepmann%rcutsq) THEN 00623 CALL siepmann_forces_v2 ( siepmann, r_last_update_pbc, cell_v, cell,& 00624 nloc_size, glob_loc_list(:,ifirst:ilast), glob_cell_v(:,ifirst:ilast),& 00625 atom_a, atom_b, f_nonbond, pv_nonbond, use_virial, siepmann%rcutsq,& 00626 particle_set, error) 00627 CALL siepmann_forces_v3( siepmann, r_last_update_pbc, cell_v,& 00628 nloc_size, glob_loc_list(:,ifirst:ilast), glob_cell_v(:,ifirst:ilast),& 00629 atom_a, atom_b, f_nonbond, pv_nonbond, use_virial, siepmann%rcutsq, & 00630 cell, particle_set, error) 00631 END IF 00632 END DO 00633 ifirst = ilast + 1 00634 IF (ipair<=npairs) junique = sort_list(1,ipair) 00635 END DO 00636 DEALLOCATE(sort_list,work_list,stat=stat) 00637 CPPostcondition(stat==0,cp_failure_level,routineP,error,failure) 00638 END IF 00639 END DO 00640 END DO Kind_Group_Loop3 00641 END DO 00642 CALL destroy_siepmann_arrays(glob_loc_list, glob_cell_v, glob_loc_list_a, error) 00643 END IF 00644 IF (use_virial) THEN 00645 pv_nonbond(1,1) = pv_nonbond(1,1) + ptens11 00646 pv_nonbond(1,2) = pv_nonbond(1,2) + ptens12 00647 pv_nonbond(1,3) = pv_nonbond(1,3) + ptens13 00648 pv_nonbond(2,1) = pv_nonbond(2,1) + ptens21 00649 pv_nonbond(2,2) = pv_nonbond(2,2) + ptens22 00650 pv_nonbond(2,3) = pv_nonbond(2,3) + ptens23 00651 pv_nonbond(3,1) = pv_nonbond(3,1) + ptens31 00652 pv_nonbond(3,2) = pv_nonbond(3,2) + ptens32 00653 pv_nonbond(3,3) = pv_nonbond(3,3) + ptens33 00654 END IF 00655 CALL timestop ( handle ) 00656 END SUBROUTINE force_nonbond_manybody 00657 00658 END MODULE manybody_potential 00659
1.7.3