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