Line data Source code
1 : !--------------------------------------------------------------------------------------------------!
2 : ! CP2K: A general program to perform molecular dynamics simulations !
3 : ! Copyright 2000-2024 CP2K developers group <https://cp2k.org> !
4 : ! !
5 : ! SPDX-License-Identifier: GPL-2.0-or-later !
6 : !--------------------------------------------------------------------------------------------------!
7 :
8 : ! **************************************************************************************************
9 : !> \brief Calculation of gCP pair potentials
10 : !> \author JGH
11 : ! **************************************************************************************************
12 : MODULE qs_gcp_method
13 : USE ai_overlap, ONLY: overlap_ab
14 : USE atomic_kind_types, ONLY: atomic_kind_type,&
15 : get_atomic_kind_set
16 : USE atprop_types, ONLY: atprop_array_init,&
17 : atprop_type
18 : USE cell_types, ONLY: cell_type
19 : USE cp_log_handling, ONLY: cp_logger_get_default_io_unit
20 : USE kinds, ONLY: dp
21 : USE message_passing, ONLY: mp_para_env_type
22 : USE particle_types, ONLY: particle_type
23 : USE physcon, ONLY: kcalmol
24 : USE qs_environment_types, ONLY: get_qs_env,&
25 : qs_environment_type
26 : USE qs_force_types, ONLY: qs_force_type
27 : USE qs_gcp_types, ONLY: qs_gcp_type
28 : USE qs_kind_types, ONLY: qs_kind_type
29 : USE qs_neighbor_list_types, ONLY: get_iterator_info,&
30 : neighbor_list_iterate,&
31 : neighbor_list_iterator_create,&
32 : neighbor_list_iterator_p_type,&
33 : neighbor_list_iterator_release,&
34 : neighbor_list_set_p_type
35 : USE virial_methods, ONLY: virial_pair_force
36 : USE virial_types, ONLY: virial_type
37 : #include "./base/base_uses.f90"
38 :
39 : IMPLICIT NONE
40 :
41 : PRIVATE
42 :
43 : CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'qs_gcp_method'
44 :
45 : PUBLIC :: calculate_gcp_pairpot
46 :
47 : ! **************************************************************************************************
48 :
49 : CONTAINS
50 :
51 : ! **************************************************************************************************
52 : !> \brief ...
53 : !> \param qs_env ...
54 : !> \param gcp_env ...
55 : !> \param energy ...
56 : !> \param calculate_forces ...
57 : !> \param ategcp ...
58 : !> \note
59 : !> \note energy_correction_type: also add gcp_env and egcp to the type
60 : !> \note
61 : ! **************************************************************************************************
62 10280 : SUBROUTINE calculate_gcp_pairpot(qs_env, gcp_env, energy, calculate_forces, ategcp)
63 :
64 : TYPE(qs_environment_type), POINTER :: qs_env
65 : TYPE(qs_gcp_type), POINTER :: gcp_env
66 : REAL(KIND=dp), INTENT(OUT) :: energy
67 : LOGICAL, INTENT(IN) :: calculate_forces
68 : REAL(KIND=dp), DIMENSION(:), OPTIONAL :: ategcp
69 :
70 : CHARACTER(LEN=*), PARAMETER :: routineN = 'calculate_gcp_pairpot'
71 :
72 : INTEGER :: atom_a, atom_b, handle, i, iatom, ikind, &
73 : jatom, jkind, mepos, natom, nkind, &
74 : nsto, unit_nr
75 10280 : INTEGER, ALLOCATABLE, DIMENSION(:) :: atom_of_kind, kind_of, ngcpat
76 : LOGICAL :: atenergy, atex, atstress, use_virial, &
77 : verbose
78 : REAL(KIND=dp) :: eama, eamb, egcp, expab, fac, fda, fdb, &
79 : gnorm, nvirta, nvirtb, rcc, sint, sqa, &
80 : sqb
81 10280 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: egcpat
82 : REAL(KIND=dp), DIMENSION(3) :: dsint, fdij, rij
83 : REAL(KIND=dp), DIMENSION(3, 3) :: dvirial
84 : REAL(KIND=dp), DIMENSION(6) :: cla, clb, rcut, zeta, zetb
85 : REAL(KIND=dp), DIMENSION(6, 6) :: sab
86 : REAL(KIND=dp), DIMENSION(6, 6, 3) :: dab
87 10280 : REAL(KIND=dp), DIMENSION(:), POINTER :: atener
88 10280 : REAL(KIND=dp), DIMENSION(:, :, :), POINTER :: atstr
89 10280 : TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
90 : TYPE(atprop_type), POINTER :: atprop
91 : TYPE(cell_type), POINTER :: cell
92 : TYPE(mp_para_env_type), POINTER :: para_env
93 : TYPE(neighbor_list_iterator_p_type), &
94 10280 : DIMENSION(:), POINTER :: nl_iterator
95 : TYPE(neighbor_list_set_p_type), DIMENSION(:), &
96 10280 : POINTER :: sab_gcp
97 10280 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
98 10280 : TYPE(qs_force_type), DIMENSION(:), POINTER :: force
99 10280 : TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
100 : TYPE(virial_type), POINTER :: virial
101 :
102 10280 : energy = 0._dp
103 10280 : IF (.NOT. gcp_env%do_gcp) RETURN
104 :
105 0 : CALL timeset(routineN, handle)
106 :
107 0 : NULLIFY (atomic_kind_set, qs_kind_set, particle_set, sab_gcp)
108 :
109 : CALL get_qs_env(qs_env=qs_env, atomic_kind_set=atomic_kind_set, qs_kind_set=qs_kind_set, &
110 0 : cell=cell, virial=virial, para_env=para_env, atprop=atprop)
111 0 : nkind = SIZE(atomic_kind_set)
112 0 : NULLIFY (particle_set)
113 0 : CALL get_qs_env(qs_env=qs_env, particle_set=particle_set)
114 0 : natom = SIZE(particle_set)
115 :
116 0 : verbose = gcp_env%verbose
117 0 : IF (verbose) THEN
118 0 : unit_nr = cp_logger_get_default_io_unit()
119 : ELSE
120 : unit_nr = -1
121 : END IF
122 : ! atomic energy and stress arrays
123 0 : atenergy = atprop%energy
124 0 : IF (atenergy) THEN
125 0 : CALL atprop_array_init(atprop%ategcp, natom)
126 0 : atener => atprop%ategcp
127 : END IF
128 0 : atstress = atprop%stress
129 0 : atstr => atprop%atstress
130 : ! external atomic energy
131 0 : atex = .FALSE.
132 0 : IF (PRESENT(ategcp)) atex = .TRUE.
133 :
134 0 : IF (unit_nr > 0) THEN
135 0 : WRITE (unit_nr, *)
136 0 : WRITE (unit_nr, *) " Pair potential geometrical counterpoise (gCP) calculation"
137 0 : WRITE (unit_nr, *)
138 0 : WRITE (unit_nr, "(T15,A,T74,F7.4)") " Gloabal Parameters: sigma = ", gcp_env%sigma, &
139 0 : " alpha = ", gcp_env%alpha, &
140 0 : " beta = ", gcp_env%beta, &
141 0 : " eta = ", gcp_env%eta
142 0 : WRITE (unit_nr, *)
143 0 : WRITE (unit_nr, "(T31,4(A5,10X))") " kind", "nvirt", "Emiss", " asto"
144 0 : DO ikind = 1, nkind
145 0 : WRITE (unit_nr, "(T31,i5,F15.1,F15.4,F15.4)") ikind, gcp_env%gcp_kind(ikind)%nbvirt, &
146 0 : gcp_env%gcp_kind(ikind)%eamiss, gcp_env%gcp_kind(ikind)%asto
147 : END DO
148 0 : WRITE (unit_nr, *)
149 : END IF
150 :
151 0 : IF (calculate_forces) THEN
152 0 : NULLIFY (force)
153 0 : CALL get_qs_env(qs_env=qs_env, force=force)
154 0 : ALLOCATE (atom_of_kind(natom), kind_of(natom))
155 0 : CALL get_atomic_kind_set(atomic_kind_set, atom_of_kind=atom_of_kind, kind_of=kind_of)
156 0 : use_virial = virial%pv_availability .AND. (.NOT. virial%pv_numer)
157 0 : IF (use_virial) dvirial = virial%pv_virial
158 : END IF
159 :
160 : ! include all integrals in the list
161 0 : rcut = 1.e6_dp
162 :
163 0 : egcp = 0.0_dp
164 0 : IF (verbose) THEN
165 0 : ALLOCATE (egcpat(natom), ngcpat(natom))
166 0 : egcpat = 0.0_dp
167 0 : ngcpat = 0
168 : END IF
169 :
170 0 : nsto = 6
171 0 : DO ikind = 1, nkind
172 0 : CPASSERT(nsto == SIZE(gcp_env%gcp_kind(jkind)%al))
173 : END DO
174 :
175 0 : sab_gcp => gcp_env%sab_gcp
176 0 : CALL neighbor_list_iterator_create(nl_iterator, sab_gcp)
177 0 : DO WHILE (neighbor_list_iterate(nl_iterator) == 0)
178 :
179 0 : CALL get_iterator_info(nl_iterator, mepos=mepos, ikind=ikind, jkind=jkind, iatom=iatom, jatom=jatom, r=rij)
180 :
181 0 : rcc = SQRT(rij(1)*rij(1) + rij(2)*rij(2) + rij(3)*rij(3))
182 0 : IF (rcc > 1.e-6_dp) THEN
183 0 : fac = 1._dp
184 0 : IF (iatom == jatom) fac = 0.5_dp
185 0 : nvirta = gcp_env%gcp_kind(ikind)%nbvirt
186 0 : nvirtb = gcp_env%gcp_kind(jkind)%nbvirt
187 0 : eama = gcp_env%gcp_kind(ikind)%eamiss
188 0 : eamb = gcp_env%gcp_kind(jkind)%eamiss
189 0 : expab = EXP(-gcp_env%alpha*rcc**gcp_env%beta)
190 0 : zeta(1:nsto) = gcp_env%gcp_kind(ikind)%al(1:nsto)
191 0 : zetb(1:nsto) = gcp_env%gcp_kind(jkind)%al(1:nsto)
192 0 : cla(1:nsto) = gcp_env%gcp_kind(ikind)%cl(1:nsto)
193 0 : clb(1:nsto) = gcp_env%gcp_kind(jkind)%cl(1:nsto)
194 0 : IF (calculate_forces) THEN
195 0 : CALL overlap_ab(0, 0, nsto, rcut, zeta, 0, 0, nsto, rcut, zetb, rij, sab, dab)
196 0 : DO i = 1, 3
197 0 : dsint(i) = SUM(cla*MATMUL(dab(:, :, i), clb))
198 : END DO
199 : ELSE
200 0 : CALL overlap_ab(0, 0, nsto, rcut, zeta, 0, 0, nsto, rcut, zetb, rij, sab)
201 : END IF
202 0 : sint = SUM(cla*MATMUL(sab, clb))
203 0 : IF (sint < 1.e-16_dp) CYCLE
204 0 : sqa = SQRT(sint*nvirta)
205 0 : sqb = SQRT(sint*nvirtb)
206 0 : IF (sqb > 1.e-12_dp) THEN
207 0 : fda = gcp_env%sigma*eama*expab/sqb
208 : ELSE
209 : fda = 0.0_dp
210 : END IF
211 0 : IF (sqa > 1.e-12_dp) THEN
212 0 : fdb = gcp_env%sigma*eamb*expab/sqa
213 : ELSE
214 : fdb = 0.0_dp
215 : END IF
216 0 : egcp = egcp + fac*(fda + fdb)
217 0 : IF (verbose) THEN
218 0 : egcpat(iatom) = egcpat(iatom) + fac*fda
219 0 : egcpat(jatom) = egcpat(jatom) + fac*fdb
220 0 : ngcpat(iatom) = ngcpat(iatom) + 1
221 0 : ngcpat(jatom) = ngcpat(jatom) + 1
222 : END IF
223 0 : IF (calculate_forces) THEN
224 0 : fdij = -fac*(fda + fdb)*(gcp_env%alpha*gcp_env%beta*rcc**(gcp_env%beta - 1.0_dp)*rij(1:3)/rcc)
225 0 : IF (sqa > 1.e-12_dp) THEN
226 0 : fdij = fdij + 0.25_dp*fac*fdb/(sqa*sqa)*dsint(1:3)
227 : END IF
228 0 : IF (sqb > 1.e-12_dp) THEN
229 0 : fdij = fdij + 0.25_dp*fac*fda/(sqb*sqb)*dsint(1:3)
230 : END IF
231 0 : atom_a = atom_of_kind(iatom)
232 0 : atom_b = atom_of_kind(jatom)
233 0 : force(ikind)%gcp(:, atom_a) = force(ikind)%gcp(:, atom_a) - fdij(:)
234 0 : force(jkind)%gcp(:, atom_b) = force(jkind)%gcp(:, atom_b) + fdij(:)
235 0 : IF (use_virial) THEN
236 0 : CALL virial_pair_force(virial%pv_virial, -1._dp, fdij, rij)
237 : END IF
238 0 : IF (atstress) THEN
239 0 : CALL virial_pair_force(atstr(:, :, iatom), -0.5_dp, fdij, rij)
240 0 : CALL virial_pair_force(atstr(:, :, jatom), -0.5_dp, fdij, rij)
241 : END IF
242 : END IF
243 0 : IF (atenergy) THEN
244 0 : atener(iatom) = atener(iatom) + fda*fac
245 0 : atener(jatom) = atener(jatom) + fdb*fac
246 : END IF
247 0 : IF (atex) THEN
248 0 : ategcp(iatom) = ategcp(iatom) + fda*fac
249 0 : ategcp(jatom) = ategcp(jatom) + fdb*fac
250 : END IF
251 : END IF
252 : END DO
253 :
254 0 : CALL neighbor_list_iterator_release(nl_iterator)
255 :
256 : ! set gCP energy
257 0 : CALL para_env%sum(egcp)
258 0 : energy = egcp
259 0 : IF (verbose) THEN
260 0 : CALL para_env%sum(egcpat)
261 0 : CALL para_env%sum(ngcpat)
262 : END IF
263 :
264 0 : IF (unit_nr > 0) THEN
265 0 : WRITE (unit_nr, "(T15,A,T61,F20.10)") " Total gCP energy [au] :", egcp
266 0 : WRITE (unit_nr, "(T15,A,T61,F20.10)") " Total gCP energy [kcal] :", egcp*kcalmol
267 0 : WRITE (unit_nr, *)
268 0 : WRITE (unit_nr, "(T19,A)") " gCP atomic energy contributions"
269 0 : WRITE (unit_nr, "(T19,A,T60,A20)") " # sites", " BSSE [kcal/mol]"
270 0 : DO i = 1, natom
271 0 : WRITE (unit_nr, "(12X,I8,10X,I8,T61,F20.10)") i, ngcpat(i), egcpat(i)*kcalmol
272 : END DO
273 : END IF
274 0 : IF (calculate_forces) THEN
275 0 : IF (unit_nr > 0) THEN
276 0 : WRITE (unit_nr, *) " gCP Forces "
277 0 : WRITE (unit_nr, *) " Atom Kind Forces "
278 : END IF
279 0 : gnorm = 0._dp
280 0 : DO iatom = 1, natom
281 0 : ikind = kind_of(iatom)
282 0 : atom_a = atom_of_kind(iatom)
283 0 : fdij(1:3) = force(ikind)%gcp(:, atom_a)
284 0 : CALL para_env%sum(fdij)
285 0 : gnorm = gnorm + SUM(ABS(fdij))
286 0 : IF (unit_nr > 0) WRITE (unit_nr, "(i5,i7,3F20.14)") iatom, ikind, fdij
287 : END DO
288 0 : IF (unit_nr > 0) THEN
289 0 : WRITE (unit_nr, *)
290 0 : WRITE (unit_nr, *) " |G| = ", gnorm
291 0 : WRITE (unit_nr, *)
292 : END IF
293 0 : IF (use_virial) THEN
294 0 : dvirial = virial%pv_virial - dvirial
295 0 : CALL para_env%sum(dvirial)
296 0 : IF (unit_nr > 0) THEN
297 0 : WRITE (unit_nr, *) " Stress Tensor (gCP)"
298 0 : WRITE (unit_nr, "(3G20.12)") dvirial
299 0 : WRITE (unit_nr, *) " Tr(P)/3 : ", (dvirial(1, 1) + dvirial(2, 2) + dvirial(3, 3))/3._dp
300 0 : WRITE (unit_nr, *)
301 : END IF
302 : END IF
303 : END IF
304 0 : IF (verbose) THEN
305 0 : DEALLOCATE (egcpat, ngcpat)
306 : END IF
307 :
308 0 : CALL timestop(handle)
309 :
310 10280 : END SUBROUTINE calculate_gcp_pairpot
311 :
312 : ! **************************************************************************************************
313 :
314 : END MODULE qs_gcp_method
|
