32 #include "./base/base_uses.f90"
36 CHARACTER(len=*),
PARAMETER,
PRIVATE :: moduleN =
'fist_efield_methods'
60 efield, use_virial, iunit, charges)
61 REAL(kind=
dp),
INTENT(OUT) :: qenergy
62 REAL(kind=
dp),
DIMENSION(:, :),
INTENT(OUT) :: qforce
63 REAL(kind=
dp),
DIMENSION(3, 3),
INTENT(OUT) :: qpv
64 TYPE(atomic_kind_type),
DIMENSION(:),
POINTER :: atomic_kind_set
65 TYPE(particle_type),
DIMENSION(:),
POINTER :: particle_set
66 TYPE(cell_type),
POINTER :: cell
68 LOGICAL,
INTENT(IN),
OPTIONAL :: use_virial
69 INTEGER,
INTENT(IN),
OPTIONAL :: iunit
70 REAL(kind=
dp),
DIMENSION(:),
OPTIONAL,
POINTER :: charges
72 COMPLEX(KIND=dp) :: zeta
73 COMPLEX(KIND=dp),
DIMENSION(3) :: ggamma
74 INTEGER :: i, ii, iparticle_kind, iw, j
75 INTEGER,
DIMENSION(:),
POINTER :: atom_list
76 LOGICAL :: use_charges, virial
77 REAL(kind=
dp) :: q, theta
78 REAL(kind=
dp),
DIMENSION(3) :: ci, dfilter, di, dipole, fieldpol, fq, &
80 TYPE(atomic_kind_type),
POINTER :: atomic_kind
87 IF (
PRESENT(charges))
THEN
88 IF (
ASSOCIATED(charges)) use_charges = .true.
91 IF (
PRESENT(iunit))
THEN
97 IF (
PRESENT(use_virial))
THEN
103 fieldpol = efield%polarisation
104 fieldpol = fieldpol/sqrt(dot_product(fieldpol, fieldpol))
105 fieldpol = -fieldpol*efield%strength
107 dfilter = efield%dfilter
110 ggamma = cmplx(1.0_dp, 0.0_dp, kind=
dp)
111 DO iparticle_kind = 1,
SIZE(atomic_kind_set)
112 atomic_kind => atomic_kind_set(iparticle_kind)
113 CALL get_atomic_kind(atomic_kind=atomic_kind, qeff=q, atom_list=atom_list)
115 DO i = 1,
SIZE(atom_list)
117 ria = particle_set(ii)%r(:)
119 IF (use_charges) q = charges(ii)
121 gvec =
twopi*cell%h_inv(j, :)
122 theta = sum(ria(:)*gvec(:))
123 zeta = cmplx(cos(theta), sin(theta), kind=
dp)**(-q)
124 ggamma(j) = ggamma(j)*zeta
130 IF (all(real(ggamma, kind=
dp) /= 0.0_dp))
THEN
131 tmp = aimag(ggamma)/real(ggamma, kind=
dp)
133 dipole = matmul(cell%hmat, ci)/
twopi
136 IF (efield%displacement)
THEN
138 di = dipole/cell%deth
140 theta = fieldpol(i) - 2._dp*
twopi*di(i)
141 qenergy = qenergy + dfilter(i)*theta**2
142 fq(i) = -dfilter(i)*theta
144 qenergy = 0.25_dp*cell%deth/
twopi*qenergy
145 DO i = 1,
SIZE(qforce, 2)
146 qforce(1:3, i) = fq(1:3)*qforce(1:3, i)
150 qenergy = sum(fieldpol*dipole)
151 DO i = 1,
SIZE(qforce, 2)
152 qforce(1:3, i) = fieldpol(1:3)*qforce(1:3, i)
157 DO iparticle_kind = 1,
SIZE(atomic_kind_set)
158 atomic_kind => atomic_kind_set(iparticle_kind)
160 DO i = 1,
SIZE(atom_list)
162 ria = particle_set(ii)%r(:)
165 qpv(j, 1:3) = qpv(j, 1:3) + qforce(j, ii)*ria(1:3)
170 IF (efield%displacement)
THEN
171 cpabort(
"Stress Tensor for constant D simulation is not working")
191 SUBROUTINE fist_dipole(fist_env, print_section, atomic_kind_set, particle_set, &
192 cell, unit_nr, charges)
193 TYPE(fist_environment_type),
POINTER :: fist_env
194 TYPE(section_vals_type),
POINTER :: print_section
195 TYPE(atomic_kind_type),
DIMENSION(:),
POINTER :: atomic_kind_set
196 TYPE(particle_type),
DIMENSION(:),
POINTER :: particle_set
197 TYPE(cell_type),
POINTER :: cell
198 INTEGER,
INTENT(IN) :: unit_nr
199 REAL(kind=
dp),
DIMENSION(:),
OPTIONAL,
POINTER :: charges
201 CHARACTER(LEN=default_string_length) :: description, dipole_type
202 COMPLEX(KIND=dp) :: dzeta, dzphase(3), zeta, zphase(3)
203 COMPLEX(KIND=dp),
DIMENSION(3) :: dggamma, ggamma
204 INTEGER :: i, iparticle_kind, j, reference
205 INTEGER,
DIMENSION(:),
POINTER :: atom_list
206 LOGICAL :: do_berry, use_charges
207 REAL(kind=
dp) :: charge_tot, ci(3), dci(3), dipole(3), dipole_deriv(3), drcc(3), dria(3), &
208 dtheta, gvec(3), q, rcc(3), ria(3), theta, tmp(3), via(3)
209 REAL(kind=
dp),
DIMENSION(:),
POINTER :: ref_point
210 TYPE(atomic_kind_type),
POINTER :: atomic_kind
211 TYPE(cp_result_type),
POINTER :: results
213 NULLIFY (atomic_kind)
215 reference =
section_get_ival(print_section, keyword_name=
"DIPOLE%REFERENCE")
217 description =
'[DIPOLE]'
220 use_charges = .false.
221 IF (
PRESENT(charges))
THEN
222 IF (
ASSOCIATED(charges)) use_charges = .true.
225 CALL get_reference_point(rcc, drcc, fist_env=fist_env, reference=reference, ref_point=ref_point)
228 dipole_deriv = 0.0_dp
231 dipole_type =
"periodic (Berry phase)"
234 IF (use_charges)
THEN
235 charge_tot = sum(charges)
237 DO i = 1,
SIZE(particle_set)
238 atomic_kind => particle_set(i)%atomic_kind
240 charge_tot = charge_tot + q
243 ria =
twopi*matmul(cell%h_inv, rcc)
244 zphase = cmplx(cos(ria), sin(ria),
dp)**charge_tot
246 dria =
twopi*matmul(cell%h_inv, drcc)
247 dzphase = charge_tot*cmplx(-sin(ria), cos(ria),
dp)**(charge_tot - 1.0_dp)*dria
249 ggamma = cmplx(1.0_dp, 0.0_dp, kind=
dp)
250 dggamma = cmplx(0.0_dp, 0.0_dp, kind=
dp)
251 DO iparticle_kind = 1,
SIZE(atomic_kind_set)
252 atomic_kind => atomic_kind_set(iparticle_kind)
253 CALL get_atomic_kind(atomic_kind=atomic_kind, qeff=q, atom_list=atom_list)
255 DO i = 1,
SIZE(atom_list)
256 ria = particle_set(atom_list(i))%r(:)
258 via = particle_set(atom_list(i))%v(:)
259 IF (use_charges) q = charges(atom_list(i))
261 gvec =
twopi*cell%h_inv(j, :)
262 theta = sum(ria(:)*gvec(:))
263 dtheta = sum(via(:)*gvec(:))
264 zeta = cmplx(cos(theta), sin(theta), kind=
dp)**(-q)
265 dzeta = -q*cmplx(-sin(theta), cos(theta), kind=
dp)**(-q - 1.0_dp)*dtheta
266 dggamma(j) = dggamma(j)*zeta + ggamma(j)*dzeta
267 ggamma(j) = ggamma(j)*zeta
271 dggamma = dggamma*zphase + ggamma*dzphase
272 ggamma = ggamma*zphase
273 IF (all(real(ggamma, kind=
dp) /= 0.0_dp))
THEN
274 tmp = aimag(ggamma)/real(ggamma, kind=
dp)
276 dci = (1.0_dp/(1.0_dp + tmp**2))* &
277 (aimag(dggamma)*real(ggamma, kind=
dp) - aimag(ggamma)*real(dggamma, kind=
dp))/(real(ggamma, kind=
dp))**2
279 dipole = matmul(cell%hmat, ci)/
twopi
280 dipole_deriv = matmul(cell%hmat, dci)/
twopi
284 CALL put_results(results, description, dipole)
286 dipole_type =
"non-periodic"
287 DO i = 1,
SIZE(particle_set)
288 atomic_kind => particle_set(i)%atomic_kind
289 ria = particle_set(i)%r(:)
292 IF (use_charges) q = charges(i)
293 dipole = dipole - q*(ria - rcc)
294 dipole_deriv(:) = dipole_deriv(:) - q*(particle_set(i)%v(:) - drcc)
298 CALL put_results(results, description, dipole)
300 IF (unit_nr > 0)
THEN
301 WRITE (unit_nr,
'(/,T2,A,T31,A50)') &
302 'MM_DIPOLE| Dipole type', adjustr(trim(dipole_type))
303 WRITE (unit_nr,
'(T2,A,T30,3(1X,F16.8))') &
304 'MM_DIPOLE| Moment [a.u.]', dipole(1:3)
305 WRITE (unit_nr,
'(T2,A,T30,3(1X,F16.8))') &
306 'MM_DIPOLE| Moment [Debye]', dipole(1:3)*
debye
307 WRITE (unit_nr,
'(T2,A,T30,3(1X,F16.8))') &
308 'MM_DIPOLE| Derivative [a.u.]', dipole_deriv(1:3)
subroutine pbc(r, r_pbc, s, s_pbc, a, b, c, alpha, beta, gamma, debug, info, pbc0, h, hinv)
...
Define the atomic kind types and their sub types.
subroutine, public get_atomic_kind(atomic_kind, fist_potential, element_symbol, name, mass, kind_number, natom, atom_list, rcov, rvdw, z, qeff, apol, cpol, mm_radius, shell, shell_active, damping)
Get attributes of an atomic kind.
Handles all functions related to the CELL.
set of type/routines to handle the storage of results in force_envs
subroutine, public cp_results_erase(results, description, nval)
erase a part of result_list
set of type/routines to handle the storage of results in force_envs
subroutine, public fist_efield_energy_force(qenergy, qforce, qpv, atomic_kind_set, particle_set, cell, efield, use_virial, iunit, charges)
...
subroutine, public fist_dipole(fist_env, print_section, atomic_kind_set, particle_set, cell, unit_nr, charges)
Evaluates the Dipole of a classical charge distribution(point-like) possibly using the berry phase fo...
subroutine, public fist_env_get(fist_env, atomic_kind_set, particle_set, ewald_pw, local_particles, local_molecules, molecule_kind_set, molecule_set, cell, cell_ref, ewald_env, fist_nonbond_env, thermo, para_env, subsys, qmmm, qmmm_env, input, shell_model, shell_model_ad, shell_particle_set, core_particle_set, multipoles, results, exclusions, efield)
Purpose: Get the FIST environment.
Defines the basic variable types.
integer, parameter, public dp
integer, parameter, public default_string_length
Definition of mathematical constants and functions.
real(kind=dp), parameter, public twopi
Calculates the moment integrals <a|r^m|b>
subroutine, public get_reference_point(rpoint, drpoint, qs_env, fist_env, reference, ref_point, ifirst, ilast)
...
Define the data structure for the particle information.
Definition of physical constants:
real(kind=dp), parameter, public debye