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