1!--------------------------------------------------------------------------------------------------! 2! DFTB+: general package for performing fast atomistic simulations ! 3! Copyright (C) 2006 - 2019 DFTB+ developers group ! 4! ! 5! See the LICENSE file for terms of usage and distribution. ! 6!--------------------------------------------------------------------------------------------------! 7 8#:include 'common.fypp' 9 10!> Routines for calculating the interaction with external charges 11module dftbp_externalcharges 12 use dftbp_assert 13 use dftbp_accuracy 14 use dftbp_blasroutines 15 use dftbp_coulomb 16 use dftbp_constants 17 use dftbp_periodic, only : getCellTranslations, foldCoordToUnitCell 18 use dftbp_environment 19 implicit none 20 21 private 22 23 public :: TExtCharge, TExtCharge_init 24 25 26 !> Private module variables 27 type TExtCharge 28 private 29 30 !> Number of point charges 31 integer :: nChrg 32 33 !> Number of atoms 34 integer :: nAtom 35 36 !> Coordinates of the point charges 37 real(dp), allocatable :: coords(:,:) 38 39 !> Charge of the point charges 40 real(dp), allocatable :: charges(:) 41 42 !> Shift vector 43 real(dp), allocatable :: invRVec(:) 44 45 !> If charges should be blured 46 logical :: tBlur 47 48 !> Blur width for the charges. 49 real(dp), allocatable :: blurWidths(:) 50 51 !> System periodic? 52 logical :: tPeriodic 53 54 !> Calculator initialised? 55 logical :: tInitialized = .false. 56 57 !> First coordinates received? 58 logical :: tUpdated = .false. 59 60 contains 61 private 62 procedure :: setCoordsCluster 63 procedure :: setCoordsPeriodic 64 procedure :: addForceDcCluster 65 procedure :: addForceDcPeriodic 66 procedure :: getElStatPotentialCluster 67 procedure :: getElStatPotentialPeriodic 68 69 !> Updates the stored coordinates for point charges 70 generic, public :: setCoordinates => setCoordsCluster, setCoordsPeriodic 71 72 !> Updates the lattice vectors 73 procedure, public :: setLatticeVectors 74 75 !> Sets the external charges 76 procedure, public :: setExternalCharges 77 78 !> Adds energy contributions per atom 79 procedure, public :: addEnergyPerAtom 80 81 !> Adds the potential from the point charges 82 procedure, public :: addShiftPerAtom 83 84 !> Adds force double counting component 85 generic, public :: addForceDc => addForceDcCluster, addForceDcPeriodic 86 87 !> Returns the electrostatic potential on a grid 88 generic, public :: getElStatPotential => getElStatPotentialCluster, getElStatPotentialPeriodic 89 90 end type TExtCharge 91 92 93contains 94 95 96 !> Initializes the calculator for external charges 97 subroutine TExtCharge_init(this, nAtom, nExtCharge, tPeriodic) 98 99 !> External charge object 100 type(TExtCharge), intent(out) :: this 101 102 !> Nr. of atoms 103 integer, intent(in) :: nAtom 104 105 !> Nr. of external charges 106 integer, intent(in) :: nExtCharge 107 108 !> Whether the system is periodic 109 logical, intent(in) :: tPeriodic 110 111 this%nAtom = nAtom 112 this%nChrg = nExtCharge 113 allocate(this%coords(3, this%nChrg)) 114 allocate(this%charges(this%nChrg)) 115 allocate(this%invRVec(nAtom)) 116 this%tPeriodic = tPeriodic 117 this%tBlur = .false. 118 119 this%tUpdated = .false. 120 this%tInitialized = .false. 121 122 end subroutine TExtCharge_init 123 124 125 !> Sets the external point charges. 126 !> 127 !> This call should be executed before setLattticeVectors() 128 !> 129 subroutine setExternalCharges(this, chargeCoords, chargeQs, blurWidths) 130 131 !> Instance 132 class(TExtCharge), intent(inout) :: this 133 134 !> Coordinates of the external charges as (3, nExtCharges) array 135 real(dp), intent(in) :: chargeCoords(:,:) 136 137 !> Charges of the point charges (sign: eletron is negative) 138 real(dp), intent(in) :: chargeQs(:) 139 140 !> Width of the Gaussians if the charges are blurred 141 real(dp), intent(in), optional :: blurWidths(:) 142 143 this%coords(:,:) = chargeCoords 144 145 ! Adapt to internal sign convention for potentials (electron has positive charge) 146 this%charges(:) = -chargeQs 147 148 if (present(blurWidths)) then 149 this%tBlur = any(blurWidths > 1.0e-7_dp) 150 else 151 this%tBlur = .false. 152 end if 153 if (this%tBlur) then 154 if (.not. allocated(this%blurWidths)) then 155 allocate(this%blurWidths(this%nChrg)) 156 end if 157 this%blurWidths(:) = blurWidths 158 end if 159 160 this%tInitialized = .true. 161 162 end subroutine setExternalCharges 163 164 165 !> Updates the module, if the lattice vectors had been changed 166 subroutine setLatticeVectors(this, latVecs, recVecs) 167 168 !> External charges structure 169 class(TExtCharge), intent(inout) :: this 170 171 !> New lattice vectors 172 real(dp), intent(in) :: latVecs(:,:) 173 174 !> New reciprocal lattice vectors. 175 real(dp), intent(in) :: recVecs(:,:) 176 177 @:ASSERT(this%tInitialized .and. this%tPeriodic) 178 179 !! Fold charges back to unit cell 180 call foldCoordToUnitCell(this%coords, latVecs, recVecs / (2.0_dp * pi)) 181 182 end subroutine setLatticeVectors 183 184 185 !> Builds the new shift vectors for new atom coordinates 186 subroutine setCoordsCluster(this, env, atomCoords) 187 188 !> External charges structure 189 class(TExtCharge), intent(inout) :: this 190 191 !> Environment settings 192 type(TEnvironment), intent(in) :: env 193 194 !> Coordinates of the atoms (not the point charges!) 195 real(dp), intent(in) :: atomCoords(:,:) 196 197 @:ASSERT(this%tInitialized) 198 @:ASSERT(.not. this%tPeriodic) 199 @:ASSERT(size(atomCoords, dim=1) == 3) 200 @:ASSERT(size(atomCoords, dim=2) >= this%nAtom) 201 202 if (this%tBlur) then 203 call sumInvR(env, this%nAtom, this%nChrg, atomCoords, this%coords, this%charges,& 204 & this%invRVec, blurWidths1=this%blurWidths) 205 else 206 call sumInvR(env, this%nAtom, this%nChrg, atomCoords, this%coords, this%charges, this%invRVec) 207 end if 208 209 this%tUpdated = .true. 210 211 end subroutine setCoordsCluster 212 213 214 !> Builds the new shift vectors for new atom coordinates 215 subroutine setCoordsPeriodic(this, env, atomCoords, rCellVec, gLat, alpha, volume) 216 217 !> External charges structure 218 class(TExtCharge), intent(inout) :: this 219 220 !> Environment settings 221 type(TEnvironment), intent(in) :: env 222 223 !> Coordinates of the atoms (not the point charges!) 224 real(dp), intent(in) :: atomCoords(:,:) 225 226 !> Real lattice points for asymmetric Ewald sum 227 real(dp), intent(in) :: rCellVec(:,:) 228 229 !> Reciprocal lattice vectors 230 real(dp), intent(in) :: gLat(:,:) 231 232 !> Ewald parameter 233 real(dp), intent(in) :: alpha 234 235 !> Cell volume 236 real(dp), intent(in) :: volume 237 238 @:ASSERT(this%tInitialized) 239 @:ASSERT(this%tPeriodic) 240 @:ASSERT(size(atomCoords, dim=1) == 3) 241 @:ASSERT(size(atomCoords, dim=2) >= this%nAtom) 242 @:ASSERT(size(gLat, dim=1) == 3) 243 244 if (this%tBlur) then 245 call sumInvR(env, this%nAtom, this%nChrg, atomCoords, this%coords, this%charges,& 246 & rCellVec, gLat, alpha, volume, this%invRVec, blurWidths1=this%blurWidths) 247 else 248 call sumInvR(env, this%nAtom, this%nChrg, atomCoords, this%coords, this%charges,& 249 & rCellVec, gLat, alpha, volume, this%invRVec) 250 end if 251 252 this%tUpdated = .true. 253 254 end subroutine setCoordsPeriodic 255 256 257 !> Adds the contribution of the external charges to the shift vector 258 subroutine addShiftPerAtom(this, shift) 259 260 !> External charges structure 261 class(TExtCharge), intent(in) :: this 262 263 !> Shift vector to add the contribution to. 264 real(dp), intent(inout) :: shift(:) 265 266 @:ASSERT(this%tInitialized .and. this%tUpdated) 267 @:ASSERT(size(shift) == this%nAtom) 268 269 shift(:) = shift + this%invRVec 270 271 end subroutine addShiftPerAtom 272 273 274 !> Adds the atomic energy contribution due to the external charges. 275 subroutine addEnergyPerAtom(this, atomCharges, energy) 276 277 !> External charges structure 278 class(TExtCharge), intent(in) :: this 279 280 !> Charge of the atoms 281 real(dp), intent(in) :: atomCharges(:) 282 283 !> Vector containing the energy per atom values. 284 real(dp), intent(inout) :: energy(:) 285 286 @:ASSERT(this%tInitialized .and. this%tUpdated) 287 @:ASSERT(size(atomCharges) == this%nAtom) 288 @:ASSERT(size(energy) == this%nAtom) 289 290 energy(:) = energy(:) + this%invRVec(:) * atomCharges(:) 291 292 end subroutine addEnergyPerAtom 293 294 295 !> Adds that part of force contribution due to the external charges, which is not contained in the 296 !> term with the shift vectors. 297 subroutine addForceDcCluster(this, env, atomForces, chrgForces, atomCoords, atomCharges) 298 299 !> External charges structure 300 class(TExtCharge), intent(in) :: this 301 302 !> Environment settings 303 type(TEnvironment), intent(in) :: env 304 305 !> Force vectors on the atoms 306 real(dp), intent(inout) :: atomForces(:,:) 307 308 !> Force vectors on the external charges 309 real(dp), intent(inout) :: chrgForces(:,:) 310 311 !> Coordinates of the atoms. 312 real(dp), intent(in) :: atomCoords(:,:) 313 314 !> Charges of the atoms. 315 real(dp), intent(in) :: atomCharges(:) 316 317 @:ASSERT(size(atomForces, dim=1) == 3) 318 @:ASSERT(size(atomForces, dim=2) == this%nAtom) 319 @:ASSERT(size(atomCoords, dim=1) == 3) 320 @:ASSERT(size(atomCoords, dim=2) == this%nAtom) 321 322 if (this%tBlur) then 323 call addInvRPrime(env, this%nAtom, this%nChrg, atomCoords, this%coords, atomCharges,& 324 & this%charges, atomForces, chrgForces, blurWidths1=this%blurWidths) 325 else 326 call addInvRPrime(env, this%nAtom, this%nChrg, atomCoords, this%coords, atomCharges,& 327 & this%charges, atomForces, chrgForces) 328 end if 329 330 end subroutine addForceDcCluster 331 332 333 !> Adds that part of force contribution due to the external charges, which is not contained in the 334 !> term with the shift vectors. 335 subroutine addForceDcPeriodic(this, env, atomForces, chrgForces, atomCoords, atomCharges,& 336 & rCellVec, gVec, alpha, vol) 337 338 !> External charges structure 339 class(TExtCharge), intent(in) :: this 340 341 !> Environment settings 342 type(TEnvironment), intent(in) :: env 343 344 !> Force vectors on the atoms 345 real(dp), intent(inout) :: atomForces(:,:) 346 347 !> Force vectors on the external charges 348 real(dp), intent(inout) :: chrgForces(:,:) 349 350 !> Coordinates of the atoms. 351 real(dp), intent(in) :: atomCoords(:,:) 352 353 !> Charges of the atoms. 354 real(dp), intent(in) :: atomCharges(:) 355 356 !> Real lattice points for asymmetric Ewald sum 357 real(dp), intent(in) :: rCellVec(:,:) 358 359 !> Reciprocal lattice vectors for the Ewald summation 360 real(dp), intent(in) :: gVec(:,:) 361 362 !> Parameter of the Ewald summation 363 real(dp), intent(in) :: alpha 364 365 !> Cell volume 366 real(dp), intent(in) :: vol 367 368 @:ASSERT(size(atomForces, dim=1) == 3) 369 @:ASSERT(size(atomForces, dim=2) == this%nAtom) 370 @:ASSERT(size(atomCoords, dim=1) == 3) 371 @:ASSERT(size(atomCoords, dim=2) >= this%nAtom) 372 373 if (this%tBlur) then 374 call addInvRPrime(env, this%nAtom, this%nChrg, atomCoords, this%coords, atomCharges,& 375 & this%charges, rCellVec, gVec, alpha, vol, atomForces, chrgForces,& 376 & blurWidths1=this%blurWidths) 377 else 378 call addInvRPrime(env, this%nAtom, this%nChrg, atomCoords, this%coords, atomCharges,& 379 & this%charges, rCellVec, gVec, alpha, vol, atomForces, chrgForces) 380 end if 381 382 end subroutine addForceDcPeriodic 383 384 385 !> Returns potential from external charges (periodic case) 386 subroutine getElStatPotentialPeriodic(this, env, locations, rCellVec, gLatPoint, alpha, volume,& 387 & V, epsSoften) 388 389 !> Instance of SCC calculation 390 class(TExtCharge), intent(in) :: this 391 392 !> Computational environment settings 393 type(TEnvironment), intent(in) :: env 394 395 !> sites to calculate potential 396 real(dp), intent(in) :: locations(:,:) 397 398 !> Real lattice points for Ewald-sum. 399 real(dp), intent(in) :: rCellVec(:,:) 400 401 !> Lattice points for reciprocal Ewald 402 real(dp), intent(in) :: gLatPoint(:,:) 403 404 !> Parameter for Ewald 405 real(dp), intent(in) :: alpha 406 407 !> Cell volume 408 real(dp), intent(in) :: volume 409 410 !> Resulting potentials 411 real(dp), intent(out) :: V(:) 412 413 !> optional potential softening 414 real(dp), optional, intent(in) :: epsSoften 415 416 @:ASSERT(all(shape(locations) == [3, size(V)])) 417 418 V(:) = 0.0_dp 419 420 if (allocated(this%blurWidths)) then 421 call sumInvR(env, size(V), size(this%charges), locations, this%coords, -this%charges,& 422 & rCellVec, gLatPoint, alpha, volume, V, this%blurWidths, epsSoften=epsSoften) 423 else 424 call sumInvR(env, size(V), size(this%charges), locations, this%coords, -this%charges,& 425 & rCellVec, gLatPoint, alpha, volume, V, epsSoften=epsSoften) 426 end if 427 428 end subroutine getElStatPotentialPeriodic 429 430 431 !> Returns potential from external charges (cluster case) 432 subroutine getElStatPotentialCluster(this, env, locations, V, epsSoften) 433 434 !> Instance of SCC calculation 435 class(TExtCharge), intent(in) :: this 436 437 !> Computational environment settings 438 type(TEnvironment), intent(in) :: env 439 440 !> sites to calculate potential 441 real(dp), intent(in) :: locations(:,:) 442 443 !> Resulting potentials 444 real(dp), intent(out) :: V(:) 445 446 !> optional potential softening 447 real(dp), optional, intent(in) :: epsSoften 448 449 @:ASSERT(all(shape(locations) == [3, size(V)])) 450 451 V(:) = 0.0_dp 452 453 if (allocated(this%blurWidths)) then 454 call sumInvR(env, size(V), size(this%charges), locations, this%coords, -this%charges, V,& 455 & this%blurWidths, epsSoften=epsSoften) 456 else 457 call sumInvR(env, size(V), size(this%charges), locations, this%coords, -this%charges, V,& 458 & epsSoften=epsSoften) 459 end if 460 461 end subroutine getElStatPotentialCluster 462 463 464end module dftbp_externalcharges 465