Line data Source code
1 : !--------------------------------------------------------------------------------------------------!
2 : ! CP2K: A general program to perform molecular dynamics simulations !
3 : ! Copyright 2000-2025 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 10538 : 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, natom, nkind, nsto, &
74 : unit_nr
75 10538 : INTEGER, ALLOCATABLE, DIMENSION(:) :: atom_of_kind, kind_of, ngcpat
76 : LOGICAL :: atenergy, atex, use_virial, verbose
77 : REAL(KIND=dp) :: eama, eamb, egcp, expab, fac, fda, fdb, &
78 : gnorm, nvirta, nvirtb, rcc, sint, sqa, &
79 : sqb
80 10538 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: egcpat
81 : REAL(KIND=dp), DIMENSION(3) :: dsint, fdij, rij
82 : REAL(KIND=dp), DIMENSION(3, 3) :: dvirial
83 : REAL(KIND=dp), DIMENSION(6) :: cla, clb, rcut, zeta, zetb
84 : REAL(KIND=dp), DIMENSION(6, 6) :: sab
85 : REAL(KIND=dp), DIMENSION(6, 6, 3) :: dab
86 10538 : REAL(KIND=dp), DIMENSION(:), POINTER :: atener
87 10538 : TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
88 : TYPE(atprop_type), POINTER :: atprop
89 : TYPE(cell_type), POINTER :: cell
90 : TYPE(mp_para_env_type), POINTER :: para_env
91 : TYPE(neighbor_list_iterator_p_type), &
92 10538 : DIMENSION(:), POINTER :: nl_iterator
93 : TYPE(neighbor_list_set_p_type), DIMENSION(:), &
94 10538 : POINTER :: sab_gcp
95 10538 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
96 10538 : TYPE(qs_force_type), DIMENSION(:), POINTER :: force
97 10538 : TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
98 : TYPE(virial_type), POINTER :: virial
99 :
100 10538 : energy = 0._dp
101 10538 : IF (.NOT. gcp_env%do_gcp) RETURN
102 :
103 4 : CALL timeset(routineN, handle)
104 :
105 4 : NULLIFY (atomic_kind_set, qs_kind_set, particle_set, sab_gcp)
106 :
107 : CALL get_qs_env(qs_env=qs_env, atomic_kind_set=atomic_kind_set, qs_kind_set=qs_kind_set, &
108 4 : cell=cell, virial=virial, para_env=para_env, atprop=atprop)
109 4 : nkind = SIZE(atomic_kind_set)
110 4 : NULLIFY (particle_set)
111 4 : CALL get_qs_env(qs_env=qs_env, particle_set=particle_set)
112 4 : natom = SIZE(particle_set)
113 :
114 4 : verbose = gcp_env%verbose
115 4 : IF (verbose) THEN
116 4 : unit_nr = cp_logger_get_default_io_unit()
117 : ELSE
118 : unit_nr = -1
119 : END IF
120 : ! atomic energy and stress arrays
121 4 : atenergy = atprop%energy
122 4 : IF (atenergy) THEN
123 0 : CALL atprop_array_init(atprop%ategcp, natom)
124 0 : atener => atprop%ategcp
125 : END IF
126 : ! external atomic energy
127 4 : atex = .FALSE.
128 4 : IF (PRESENT(ategcp)) atex = .TRUE.
129 :
130 4 : IF (unit_nr > 0) THEN
131 2 : WRITE (unit_nr, *)
132 2 : WRITE (unit_nr, *) " Pair potential geometrical counterpoise (gCP) calculation"
133 2 : WRITE (unit_nr, *)
134 2 : WRITE (unit_nr, "(T15,A,T74,F7.4)") " Gloabal Parameters: sigma = ", gcp_env%sigma, &
135 2 : " alpha = ", gcp_env%alpha, &
136 2 : " beta = ", gcp_env%beta, &
137 4 : " eta = ", gcp_env%eta
138 2 : WRITE (unit_nr, *)
139 2 : WRITE (unit_nr, "(T31,4(A5,10X))") " kind", "nvirt", "Emiss", " asto"
140 5 : DO ikind = 1, nkind
141 3 : WRITE (unit_nr, "(T31,i5,F15.1,F15.4,F15.4)") ikind, gcp_env%gcp_kind(ikind)%nbvirt, &
142 8 : gcp_env%gcp_kind(ikind)%eamiss, gcp_env%gcp_kind(ikind)%asto
143 : END DO
144 2 : WRITE (unit_nr, *)
145 : END IF
146 :
147 4 : IF (calculate_forces) THEN
148 2 : NULLIFY (force)
149 2 : CALL get_qs_env(qs_env=qs_env, force=force)
150 2 : CALL get_atomic_kind_set(atomic_kind_set, atom_of_kind=atom_of_kind, kind_of=kind_of)
151 2 : use_virial = virial%pv_availability .AND. (.NOT. virial%pv_numer)
152 0 : IF (use_virial) dvirial = virial%pv_virial
153 : END IF
154 :
155 : ! include all integrals in the list
156 28 : rcut = 1.e6_dp
157 :
158 4 : egcp = 0.0_dp
159 4 : IF (verbose) THEN
160 20 : ALLOCATE (egcpat(natom), ngcpat(natom))
161 14 : egcpat = 0.0_dp
162 14 : ngcpat = 0
163 : END IF
164 :
165 4 : nsto = 6
166 10 : DO ikind = 1, nkind
167 4 : CPASSERT(nsto == SIZE(gcp_env%gcp_kind(jkind)%al))
168 : END DO
169 :
170 4 : sab_gcp => gcp_env%sab_gcp
171 4 : CALL neighbor_list_iterator_create(nl_iterator, sab_gcp)
172 13 : DO WHILE (neighbor_list_iterate(nl_iterator) == 0)
173 :
174 9 : CALL get_iterator_info(nl_iterator, ikind=ikind, jkind=jkind, iatom=iatom, jatom=jatom, r=rij)
175 :
176 9 : rcc = SQRT(rij(1)*rij(1) + rij(2)*rij(2) + rij(3)*rij(3))
177 13 : IF (rcc > 1.e-6_dp) THEN
178 4 : fac = 1._dp
179 4 : IF (iatom == jatom) fac = 0.5_dp
180 4 : nvirta = gcp_env%gcp_kind(ikind)%nbvirt
181 4 : nvirtb = gcp_env%gcp_kind(jkind)%nbvirt
182 4 : eama = gcp_env%gcp_kind(ikind)%eamiss
183 4 : eamb = gcp_env%gcp_kind(jkind)%eamiss
184 4 : expab = EXP(-gcp_env%alpha*rcc**gcp_env%beta)
185 28 : zeta(1:nsto) = gcp_env%gcp_kind(ikind)%al(1:nsto)
186 28 : zetb(1:nsto) = gcp_env%gcp_kind(jkind)%al(1:nsto)
187 28 : cla(1:nsto) = gcp_env%gcp_kind(ikind)%cl(1:nsto)
188 28 : clb(1:nsto) = gcp_env%gcp_kind(jkind)%cl(1:nsto)
189 4 : IF (calculate_forces) THEN
190 3 : CALL overlap_ab(0, 0, nsto, rcut, zeta, 0, 0, nsto, rcut, zetb, rij, sab, dab)
191 12 : DO i = 1, 3
192 444 : dsint(i) = SUM(cla*MATMUL(dab(:, :, i), clb))
193 : END DO
194 : ELSE
195 1 : CALL overlap_ab(0, 0, nsto, rcut, zeta, 0, 0, nsto, rcut, zetb, rij, sab)
196 : END IF
197 196 : sint = SUM(cla*MATMUL(sab, clb))
198 4 : IF (sint < 1.e-16_dp) CYCLE
199 4 : sqa = SQRT(sint*nvirta)
200 4 : sqb = SQRT(sint*nvirtb)
201 4 : IF (sqb > 1.e-12_dp) THEN
202 4 : fda = gcp_env%sigma*eama*expab/sqb
203 : ELSE
204 : fda = 0.0_dp
205 : END IF
206 4 : IF (sqa > 1.e-12_dp) THEN
207 4 : fdb = gcp_env%sigma*eamb*expab/sqa
208 : ELSE
209 : fdb = 0.0_dp
210 : END IF
211 4 : egcp = egcp + fac*(fda + fdb)
212 4 : IF (verbose) THEN
213 4 : egcpat(iatom) = egcpat(iatom) + fac*fda
214 4 : egcpat(jatom) = egcpat(jatom) + fac*fdb
215 4 : ngcpat(iatom) = ngcpat(iatom) + 1
216 4 : ngcpat(jatom) = ngcpat(jatom) + 1
217 : END IF
218 4 : IF (calculate_forces) THEN
219 12 : fdij = -fac*(fda + fdb)*(gcp_env%alpha*gcp_env%beta*rcc**(gcp_env%beta - 1.0_dp)*rij(1:3)/rcc)
220 3 : IF (sqa > 1.e-12_dp) THEN
221 12 : fdij = fdij + 0.25_dp*fac*fdb/(sqa*sqa)*dsint(1:3)
222 : END IF
223 3 : IF (sqb > 1.e-12_dp) THEN
224 12 : fdij = fdij + 0.25_dp*fac*fda/(sqb*sqb)*dsint(1:3)
225 : END IF
226 3 : atom_a = atom_of_kind(iatom)
227 3 : atom_b = atom_of_kind(jatom)
228 12 : force(ikind)%gcp(:, atom_a) = force(ikind)%gcp(:, atom_a) - fdij(:)
229 12 : force(jkind)%gcp(:, atom_b) = force(jkind)%gcp(:, atom_b) + fdij(:)
230 3 : IF (use_virial) THEN
231 0 : CALL virial_pair_force(virial%pv_virial, -1._dp, fdij, rij)
232 : END IF
233 : END IF
234 4 : IF (atenergy) THEN
235 0 : atener(iatom) = atener(iatom) + fda*fac
236 0 : atener(jatom) = atener(jatom) + fdb*fac
237 : END IF
238 4 : IF (atex) THEN
239 0 : ategcp(iatom) = ategcp(iatom) + fda*fac
240 0 : ategcp(jatom) = ategcp(jatom) + fdb*fac
241 : END IF
242 : END IF
243 : END DO
244 :
245 4 : CALL neighbor_list_iterator_release(nl_iterator)
246 :
247 : ! set gCP energy
248 4 : CALL para_env%sum(egcp)
249 4 : energy = egcp
250 4 : IF (verbose) THEN
251 4 : CALL para_env%sum(egcpat)
252 4 : CALL para_env%sum(ngcpat)
253 : END IF
254 :
255 4 : IF (unit_nr > 0) THEN
256 2 : WRITE (unit_nr, "(T15,A,T61,F20.10)") " Total gCP energy [au] :", egcp
257 2 : WRITE (unit_nr, "(T15,A,T61,F20.10)") " Total gCP energy [kcal] :", egcp*kcalmol
258 2 : WRITE (unit_nr, *)
259 2 : WRITE (unit_nr, "(T19,A)") " gCP atomic energy contributions"
260 2 : WRITE (unit_nr, "(T19,A,T60,A20)") " # sites", " BSSE [kcal/mol]"
261 7 : DO i = 1, natom
262 7 : WRITE (unit_nr, "(12X,I8,10X,I8,T61,F20.10)") i, ngcpat(i), egcpat(i)*kcalmol
263 : END DO
264 : END IF
265 4 : IF (calculate_forces) THEN
266 2 : IF (unit_nr > 0) THEN
267 1 : WRITE (unit_nr, *) " gCP Forces "
268 1 : WRITE (unit_nr, *) " Atom Kind Forces "
269 : END IF
270 2 : gnorm = 0._dp
271 8 : DO iatom = 1, natom
272 6 : ikind = kind_of(iatom)
273 6 : atom_a = atom_of_kind(iatom)
274 24 : fdij(1:3) = force(ikind)%gcp(:, atom_a)
275 6 : CALL para_env%sum(fdij)
276 24 : gnorm = gnorm + SUM(ABS(fdij))
277 8 : IF (unit_nr > 0) WRITE (unit_nr, "(i5,i7,3F20.14)") iatom, ikind, fdij
278 : END DO
279 2 : IF (unit_nr > 0) THEN
280 1 : WRITE (unit_nr, *)
281 1 : WRITE (unit_nr, *) " |G| = ", gnorm
282 1 : WRITE (unit_nr, *)
283 : END IF
284 2 : IF (use_virial) THEN
285 0 : dvirial = virial%pv_virial - dvirial
286 0 : CALL para_env%sum(dvirial)
287 0 : IF (unit_nr > 0) THEN
288 0 : WRITE (unit_nr, *) " Stress Tensor (gCP)"
289 0 : WRITE (unit_nr, "(3G20.12)") dvirial
290 0 : WRITE (unit_nr, *) " Tr(P)/3 : ", (dvirial(1, 1) + dvirial(2, 2) + dvirial(3, 3))/3._dp
291 0 : WRITE (unit_nr, *)
292 : END IF
293 : END IF
294 : END IF
295 4 : IF (verbose) THEN
296 4 : DEALLOCATE (egcpat, ngcpat)
297 : END IF
298 :
299 4 : CALL timestop(handle)
300 :
301 10542 : END SUBROUTINE calculate_gcp_pairpot
302 :
303 : ! **************************************************************************************************
304 :
305 : END MODULE qs_gcp_method
|
