1!--------------------------------------------------------------------------------------------------! 2! CP2K: A general program to perform molecular dynamics simulations ! 3! Copyright (C) 2000 - 2019 CP2K developers group ! 4!--------------------------------------------------------------------------------------------------! 5 6! ************************************************************************************************** 7!> \par History 8!> EAM potential 9!> \author CJM, I-Feng W. Kuo, Teodoro Laino 10! ************************************************************************************************** 11MODULE manybody_eam 12 13 USE bibliography, ONLY: Foiles1986,& 14 cite_reference 15 USE cell_types, ONLY: cell_type 16 USE cp_para_types, ONLY: cp_para_env_type 17 USE fist_neighbor_list_types, ONLY: fist_neighbor_type,& 18 neighbor_kind_pairs_type 19 USE fist_nonbond_env_types, ONLY: eam_type,& 20 fist_nonbond_env_get,& 21 fist_nonbond_env_set,& 22 fist_nonbond_env_type,& 23 pos_type 24 USE kinds, ONLY: dp 25 USE mathlib, ONLY: matvec_3x3 26 USE message_passing, ONLY: mp_sum 27 USE pair_potential_types, ONLY: ea_type,& 28 eam_pot_type,& 29 pair_potential_pp_type 30 USE particle_types, ONLY: particle_type 31#include "./base/base_uses.f90" 32 33 IMPLICIT NONE 34 35 PRIVATE 36 PUBLIC :: get_force_eam, density_nonbond 37 CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'manybody_eam' 38 39CONTAINS 40 41! ************************************************************************************************** 42!> \brief ... 43!> \param fist_nonbond_env ... 44!> \param particle_set ... 45!> \param cell ... 46!> \param para_env ... 47!> \author CJM 48! ************************************************************************************************** 49 SUBROUTINE density_nonbond(fist_nonbond_env, particle_set, cell, para_env) 50 TYPE(fist_nonbond_env_type), POINTER :: fist_nonbond_env 51 TYPE(particle_type), DIMENSION(:), INTENT(INOUT) :: particle_set 52 TYPE(cell_type), POINTER :: cell 53 TYPE(cp_para_env_type), POINTER :: para_env 54 55 CHARACTER(LEN=*), PARAMETER :: routineN = 'density_nonbond', & 56 routineP = moduleN//':'//routineN 57 58 INTEGER :: atom_a, atom_b, handle, i, i_a, i_b, & 59 iend, igrp, ikind, ilist, ipair, & 60 iparticle, istart, jkind, kind_a, & 61 kind_b, nkinds, nparticle 62 INTEGER, ALLOCATABLE, DIMENSION(:, :) :: eam_kinds_index 63 LOGICAL :: do_eam 64 REAL(KIND=dp) :: fac, rab2, rab2_max 65 REAL(KIND=dp), DIMENSION(3) :: cell_v, rab 66 REAL(KIND=dp), DIMENSION(:), POINTER :: rho 67 TYPE(eam_pot_type), POINTER :: eam_a, eam_b 68 TYPE(eam_type), DIMENSION(:), POINTER :: eam_data 69 TYPE(fist_neighbor_type), POINTER :: nonbonded 70 TYPE(neighbor_kind_pairs_type), POINTER :: neighbor_kind_pair 71 TYPE(pair_potential_pp_type), POINTER :: potparm 72 TYPE(pos_type), DIMENSION(:), POINTER :: r_last_update, r_last_update_pbc 73 74 CALL timeset(routineN, handle) 75 do_eam = .FALSE. 76 CALL fist_nonbond_env_get(fist_nonbond_env, nonbonded=nonbonded, & 77 potparm=potparm, r_last_update=r_last_update, & 78 r_last_update_pbc=r_last_update_pbc, eam_data=eam_data) 79 nkinds = SIZE(potparm%pot, 1) 80 ALLOCATE (eam_kinds_index(nkinds, nkinds)) 81 eam_kinds_index = -1 82 DO ikind = 1, nkinds 83 DO jkind = ikind, nkinds 84 DO i = 1, SIZE(potparm%pot(ikind, jkind)%pot%type) 85 IF (potparm%pot(ikind, jkind)%pot%type(i) == ea_type) THEN 86 ! At the moment we allow only 1 EAM per each kinds pair.. 87 CPASSERT(eam_kinds_index(ikind, jkind) == -1) 88 CPASSERT(eam_kinds_index(jkind, ikind) == -1) 89 eam_kinds_index(ikind, jkind) = i 90 eam_kinds_index(jkind, ikind) = i 91 do_eam = .TRUE. 92 END IF 93 END DO 94 END DO 95 END DO 96 97 nparticle = SIZE(particle_set) 98 99 IF (do_eam) THEN 100 IF (.NOT. ASSOCIATED(eam_data)) THEN 101 ALLOCATE (eam_data(nparticle)) 102 CALL fist_nonbond_env_set(fist_nonbond_env, eam_data=eam_data) 103 ENDIF 104 DO i = 1, nparticle 105 eam_data(i)%rho = 0.0_dp 106 eam_data(i)%f_embed = 0.0_dp 107 ENDDO 108 ENDIF 109 110 ! Only if EAM potential are present 111 IF (do_eam) THEN 112 ! Add citation 113 CALL cite_reference(Foiles1986) 114 NULLIFY (eam_a, eam_b) 115 ALLOCATE (rho(nparticle)) 116 rho = 0._dp 117 ! Starting the force loop 118 DO ilist = 1, nonbonded%nlists 119 neighbor_kind_pair => nonbonded%neighbor_kind_pairs(ilist) 120 IF (neighbor_kind_pair%npairs == 0) CYCLE 121 Kind_Group_Loop: DO igrp = 1, neighbor_kind_pair%ngrp_kind 122 istart = neighbor_kind_pair%grp_kind_start(igrp) 123 iend = neighbor_kind_pair%grp_kind_end(igrp) 124 ikind = neighbor_kind_pair%ij_kind(1, igrp) 125 jkind = neighbor_kind_pair%ij_kind(2, igrp) 126 127 i = eam_kinds_index(ikind, jkind) 128 IF (i == -1) CYCLE 129 rab2_max = potparm%pot(ikind, jkind)%pot%rcutsq 130 CALL matvec_3x3(cell_v, cell%hmat, REAL(neighbor_kind_pair%cell_vector, KIND=dp)) 131 DO ipair = istart, iend 132 atom_a = neighbor_kind_pair%list(1, ipair) 133 atom_b = neighbor_kind_pair%list(2, ipair) 134 fac = 1.0_dp 135 IF (atom_a == atom_b) fac = 0.5_dp 136 rab = r_last_update_pbc(atom_b)%r - r_last_update_pbc(atom_a)%r 137 rab = rab + cell_v 138 rab2 = rab(1)*rab(1) + rab(2)*rab(2) + rab(3)*rab(3) 139 IF (rab2 <= rab2_max) THEN 140 kind_a = particle_set(atom_a)%atomic_kind%kind_number 141 kind_b = particle_set(atom_b)%atomic_kind%kind_number 142 i_a = eam_kinds_index(kind_a, kind_a) 143 i_b = eam_kinds_index(kind_b, kind_b) 144 eam_a => potparm%pot(kind_a, kind_a)%pot%set(i_a)%eam 145 eam_b => potparm%pot(kind_b, kind_b)%pot%set(i_b)%eam 146 CALL get_rho_eam(eam_a, eam_b, rab2, atom_a, atom_b, rho, fac) 147 END IF 148 END DO 149 END DO Kind_Group_Loop 150 END DO 151 CALL mp_sum(rho, para_env%group) 152 DO iparticle = 1, nparticle 153 eam_data(iparticle)%rho = rho(iparticle) 154 END DO 155 156 DEALLOCATE (rho) 157 END IF 158 DEALLOCATE (eam_kinds_index) 159 CALL timestop(handle) 160 161 END SUBROUTINE density_nonbond 162 163! ************************************************************************************************** 164!> \brief ... 165!> \param eam_a ... 166!> \param eam_b ... 167!> \param rab2 ... 168!> \param atom_a ... 169!> \param atom_b ... 170!> \param rho ... 171!> \param fac ... 172!> \author CJM 173! ************************************************************************************************** 174 SUBROUTINE get_rho_eam(eam_a, eam_b, rab2, atom_a, atom_b, rho, fac) 175 TYPE(eam_pot_type), POINTER :: eam_a, eam_b 176 REAL(dp), INTENT(IN) :: rab2 177 INTEGER, INTENT(IN) :: atom_a, atom_b 178 REAL(dp), INTENT(INOUT) :: rho(:) 179 REAL(dp), INTENT(IN) :: fac 180 181 CHARACTER(LEN=*), PARAMETER :: routineN = 'get_rho_eam', routineP = moduleN//':'//routineN 182 183 INTEGER :: index 184 REAL(dp) :: qq, rab, rhoi, rhoj 185 186 rab = SQRT(rab2) 187 188 ! Particle A 189 index = INT(rab/eam_b%drar) + 1 190 IF (index > eam_b%npoints) THEN 191 index = eam_b%npoints 192 ELSEIF (index < 1) THEN 193 index = 1 194 ENDIF 195 qq = rab - eam_b%rval(index) 196 rhoi = eam_b%rho(index) + qq*eam_b%rhop(index) 197 198 ! Particle B 199 index = INT(rab/eam_a%drar) + 1 200 IF (index > eam_a%npoints) THEN 201 index = eam_a%npoints 202 ELSEIF (index < 1) THEN 203 index = 1 204 ENDIF 205 qq = rab - eam_a%rval(index) 206 rhoj = eam_a%rho(index) + qq*eam_a%rhop(index) 207 208 rho(atom_a) = rho(atom_a) + rhoi*fac 209 rho(atom_b) = rho(atom_b) + rhoj*fac 210 END SUBROUTINE get_rho_eam 211 212! ************************************************************************************************** 213!> \brief ... 214!> \param rab2 ... 215!> \param eam_a ... 216!> \param eam_b ... 217!> \param eam_data ... 218!> \param atom_a ... 219!> \param atom_b ... 220!> \param f_eam ... 221!> \author CJM 222! ************************************************************************************************** 223 SUBROUTINE get_force_eam(rab2, eam_a, eam_b, eam_data, atom_a, atom_b, f_eam) 224 REAL(dp), INTENT(IN) :: rab2 225 TYPE(eam_pot_type), POINTER :: eam_a, eam_b 226 TYPE(eam_type), INTENT(IN) :: eam_data(:) 227 INTEGER, INTENT(IN) :: atom_a, atom_b 228 REAL(dp), INTENT(OUT) :: f_eam 229 230 CHARACTER(LEN=*), PARAMETER :: routineN = 'get_force_eam', routineP = moduleN//':'//routineN 231 232 INTEGER :: index 233 REAL(KIND=dp) :: denspi, denspj, fcp, qq, rab 234 235 rab = SQRT(rab2) 236 237 ! Particle A 238 index = INT(rab/eam_a%drar) + 1 239 IF (index > eam_a%npoints) THEN 240 index = eam_a%npoints 241 ELSEIF (index < 1) THEN 242 index = 1 243 ENDIF 244 qq = rab - eam_a%rval(index) 245 IF (index == eam_a%npoints) THEN 246 denspi = eam_a%rhop(index) + qq*(eam_a%rhop(index) - eam_a%rhop(index - 1))/eam_a%drar 247 ELSE 248 denspi = eam_a%rhop(index) + qq*(eam_a%rhop(index + 1) - eam_a%rhop(index))/eam_a%drar 249 END IF 250 251 ! Particle B 252 index = INT(rab/eam_b%drar) + 1 253 IF (index > eam_b%npoints) THEN 254 index = eam_b%npoints 255 ELSEIF (index < 1) THEN 256 index = 1 257 ENDIF 258 qq = rab - eam_b%rval(index) 259 IF (index == eam_b%npoints) THEN 260 denspj = eam_b%rhop(index) + qq*(eam_b%rhop(index) - eam_b%rhop(index - 1))/eam_b%drar 261 ELSE 262 denspj = eam_b%rhop(index) + qq*(eam_b%rhop(index + 1) - eam_b%rhop(index))/eam_b%drar 263 END IF 264 265 fcp = denspj*eam_data(atom_a)%f_embed + denspi*eam_data(atom_b)%f_embed 266 f_eam = fcp/rab 267 END SUBROUTINE get_force_eam 268 269END MODULE manybody_eam 270 271