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 QMMM Coulomb contributions in TB
10 : !> \author JGH
11 : ! **************************************************************************************************
12 : MODULE qmmm_tb_coulomb
13 : USE atomic_kind_types, ONLY: atomic_kind_type,&
14 : get_atomic_kind_set
15 : USE atprop_types, ONLY: atprop_type
16 : USE cell_types, ONLY: cell_type,&
17 : get_cell
18 : USE cp_control_types, ONLY: dft_control_type
19 : USE cp_dbcsr_api, ONLY: dbcsr_add,&
20 : dbcsr_get_block_p,&
21 : dbcsr_iterator_blocks_left,&
22 : dbcsr_iterator_next_block,&
23 : dbcsr_iterator_start,&
24 : dbcsr_iterator_stop,&
25 : dbcsr_iterator_type,&
26 : dbcsr_p_type
27 : USE distribution_1d_types, ONLY: distribution_1d_type
28 : USE ewald_environment_types, ONLY: ewald_env_get,&
29 : ewald_environment_type
30 : USE ewald_methods_tb, ONLY: tb_spme_evaluate
31 : USE ewald_pw_types, ONLY: ewald_pw_type
32 : USE kinds, ONLY: dp
33 : USE mathconstants, ONLY: oorootpi,&
34 : pi
35 : USE message_passing, ONLY: mp_para_env_type
36 : USE particle_types, ONLY: particle_type
37 : USE pw_poisson_types, ONLY: do_ewald_ewald,&
38 : do_ewald_none,&
39 : do_ewald_pme,&
40 : do_ewald_spme
41 : USE qmmm_types_low, ONLY: qmmm_env_qm_type
42 : USE qs_energy_types, ONLY: qs_energy_type
43 : USE qs_environment_types, ONLY: get_qs_env,&
44 : qs_environment_type
45 : USE qs_force_types, ONLY: qs_force_type
46 : USE qs_rho_types, ONLY: qs_rho_get,&
47 : qs_rho_type
48 : USE virial_types, ONLY: virial_type
49 : #include "./base/base_uses.f90"
50 :
51 : IMPLICIT NONE
52 :
53 : PRIVATE
54 :
55 : CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'qmmm_tb_coulomb'
56 :
57 : PUBLIC :: build_tb_coulomb_qmqm
58 :
59 : CONTAINS
60 :
61 : ! **************************************************************************************************
62 : !> \brief ...
63 : !> \param qs_env ...
64 : !> \param ks_matrix ...
65 : !> \param rho ...
66 : !> \param mcharge ...
67 : !> \param energy ...
68 : !> \param calculate_forces ...
69 : !> \param just_energy ...
70 : ! **************************************************************************************************
71 2488 : SUBROUTINE build_tb_coulomb_qmqm(qs_env, ks_matrix, rho, mcharge, energy, &
72 : calculate_forces, just_energy)
73 :
74 : TYPE(qs_environment_type), INTENT(IN) :: qs_env
75 : TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: ks_matrix
76 : TYPE(qs_rho_type), POINTER :: rho
77 : REAL(dp), DIMENSION(:) :: mcharge
78 : TYPE(qs_energy_type), POINTER :: energy
79 : LOGICAL, INTENT(in) :: calculate_forces, just_energy
80 :
81 : CHARACTER(len=*), PARAMETER :: routineN = 'build_tb_coulomb_qmqm'
82 :
83 : INTEGER :: atom_i, atom_j, blk, ewald_type, handle, &
84 : i, ia, iatom, ikind, jatom, jkind, &
85 : natom, nmat
86 2488 : INTEGER, ALLOCATABLE, DIMENSION(:) :: atom_of_kind, kind_of
87 : INTEGER, DIMENSION(3) :: periodic
88 : LOGICAL :: found, use_virial
89 : REAL(KIND=dp) :: alpha, deth, dfr, dr, fi, fr, gmij
90 : REAL(KIND=dp), DIMENSION(3) :: fij, rij
91 2488 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: dsblock, gmcharge, ksblock, ksblock_2, &
92 2488 : pblock, sblock
93 2488 : TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
94 : TYPE(atprop_type), POINTER :: atprop
95 : TYPE(cell_type), POINTER :: cell, mm_cell
96 : TYPE(dbcsr_iterator_type) :: iter
97 2488 : TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_p, matrix_s
98 : TYPE(dft_control_type), POINTER :: dft_control
99 : TYPE(distribution_1d_type), POINTER :: local_particles
100 : TYPE(ewald_environment_type), POINTER :: ewald_env
101 : TYPE(ewald_pw_type), POINTER :: ewald_pw
102 : TYPE(mp_para_env_type), POINTER :: para_env
103 2488 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
104 : TYPE(qmmm_env_qm_type), POINTER :: qmmm_env_qm
105 2488 : TYPE(qs_force_type), DIMENSION(:), POINTER :: force
106 : TYPE(virial_type), POINTER :: virial
107 :
108 2488 : CALL timeset(routineN, handle)
109 :
110 2488 : NULLIFY (matrix_p, matrix_s, virial, atprop, dft_control)
111 :
112 2488 : use_virial = .FALSE.
113 :
114 2488 : IF (calculate_forces) THEN
115 : nmat = 4
116 : ELSE
117 2480 : nmat = 1
118 : END IF
119 :
120 2488 : natom = SIZE(mcharge)
121 9952 : ALLOCATE (gmcharge(natom, nmat))
122 12536 : gmcharge = 0._dp
123 :
124 : CALL get_qs_env(qs_env, &
125 : particle_set=particle_set, &
126 : cell=cell, &
127 : virial=virial, &
128 : atprop=atprop, &
129 2488 : dft_control=dft_control)
130 :
131 2488 : IF (calculate_forces) THEN
132 16 : use_virial = virial%pv_availability .AND. (.NOT. virial%pv_numer)
133 : END IF
134 :
135 : ! Qm-QM long range correction for QMMM calculations
136 : ! no atomic energy evaluation
137 2488 : CPASSERT(.NOT. atprop%energy)
138 : ! no stress tensor possible for QMMM
139 2488 : CPASSERT(.NOT. use_virial)
140 2488 : qmmm_env_qm => qs_env%qmmm_env_qm
141 2488 : ewald_env => qmmm_env_qm%ewald_env
142 2488 : ewald_pw => qmmm_env_qm%ewald_pw
143 2488 : CALL get_qs_env(qs_env=qs_env, super_cell=mm_cell)
144 2488 : CALL get_cell(cell=mm_cell, periodic=periodic, deth=deth)
145 2488 : CALL ewald_env_get(ewald_env, alpha=alpha, ewald_type=ewald_type)
146 12536 : gmcharge = 0.0_dp
147 : ! direct sum for overlap and local correction
148 : CALL get_qs_env(qs_env=qs_env, &
149 : atomic_kind_set=atomic_kind_set, &
150 : local_particles=local_particles, &
151 2488 : force=force, para_env=para_env)
152 7464 : DO ikind = 1, SIZE(local_particles%n_el)
153 11196 : DO ia = 1, local_particles%n_el(ikind)
154 3732 : iatom = local_particles%list(ikind)%array(ia)
155 12440 : DO jatom = 1, iatom - 1
156 14928 : rij = particle_set(iatom)%r - particle_set(jatom)%r
157 : ! no pbc(rij,mm_cell) at this point, we assume that QM particles are
158 : ! inside QM box and QM box << MM box
159 14928 : dr = SQRT(SUM(rij(:)**2))
160 : ! local (unit cell) correction 1/R
161 3732 : gmcharge(iatom, 1) = gmcharge(iatom, 1) - mcharge(jatom)/dr
162 3732 : gmcharge(jatom, 1) = gmcharge(jatom, 1) - mcharge(iatom)/dr
163 3768 : DO i = 2, nmat
164 36 : gmcharge(iatom, i) = gmcharge(iatom, i) - rij(i - 1)*mcharge(jatom)/dr**3
165 3768 : gmcharge(jatom, i) = gmcharge(jatom, i) + rij(i - 1)*mcharge(iatom)/dr**3
166 : END DO
167 : ! overlap correction
168 3732 : fr = erfc(alpha*dr)/dr
169 3732 : gmcharge(iatom, 1) = gmcharge(iatom, 1) + mcharge(jatom)*fr
170 3732 : gmcharge(jatom, 1) = gmcharge(jatom, 1) + mcharge(iatom)*fr
171 7464 : IF (nmat > 1) THEN
172 12 : dfr = -2._dp*alpha*EXP(-alpha*alpha*dr*dr)*oorootpi/dr - fr/dr
173 12 : dfr = -dfr/dr
174 48 : DO i = 2, nmat
175 36 : gmcharge(iatom, i) = gmcharge(iatom, i) - rij(i - 1)*mcharge(jatom)*dfr
176 48 : gmcharge(jatom, i) = gmcharge(jatom, i) + rij(i - 1)*mcharge(iatom)*dfr
177 : END DO
178 : END IF
179 : END DO
180 : END DO
181 : END DO
182 :
183 0 : SELECT CASE (ewald_type)
184 : CASE DEFAULT
185 0 : CPABORT("Invalid Ewald type")
186 : CASE (do_ewald_none)
187 0 : CPABORT("Not allowed with DFTB")
188 : CASE (do_ewald_ewald)
189 0 : CPABORT("Standard Ewald not implemented in DFTB")
190 : CASE (do_ewald_pme)
191 0 : CPABORT("PME not implemented in DFTB")
192 : CASE (do_ewald_spme)
193 : CALL tb_spme_evaluate(ewald_env, ewald_pw, particle_set, mm_cell, &
194 2488 : gmcharge, mcharge, calculate_forces, virial, use_virial)
195 : END SELECT
196 : !
197 17416 : CALL para_env%sum(gmcharge(:, 1))
198 : !
199 : ! add self charge interaction and background charge contribution
200 9952 : gmcharge(:, 1) = gmcharge(:, 1) - 2._dp*alpha*oorootpi*mcharge(:)
201 2488 : IF (ANY(periodic(:) == 1)) THEN
202 9952 : gmcharge(:, 1) = gmcharge(:, 1) - pi/alpha**2/deth
203 : END IF
204 : !
205 9952 : energy%qmmm_el = energy%qmmm_el + 0.5_dp*SUM(mcharge(:)*gmcharge(:, 1))
206 : !
207 2488 : IF (calculate_forces) THEN
208 : CALL get_atomic_kind_set(atomic_kind_set=atomic_kind_set, &
209 : kind_of=kind_of, &
210 8 : atom_of_kind=atom_of_kind)
211 : END IF
212 : !
213 2488 : IF (.NOT. just_energy) THEN
214 2488 : CALL get_qs_env(qs_env=qs_env, matrix_s=matrix_s)
215 2488 : CALL qs_rho_get(rho, rho_ao=matrix_p)
216 :
217 2488 : IF (calculate_forces .AND. SIZE(matrix_p) == 2) THEN
218 : CALL dbcsr_add(matrix_p(1)%matrix, matrix_p(2)%matrix, &
219 0 : alpha_scalar=1.0_dp, beta_scalar=1.0_dp)
220 : END IF
221 :
222 2488 : CALL dbcsr_iterator_start(iter, ks_matrix(1, 1)%matrix)
223 9952 : DO WHILE (dbcsr_iterator_blocks_left(iter))
224 7464 : CALL dbcsr_iterator_next_block(iter, iatom, jatom, ksblock, blk)
225 7464 : NULLIFY (sblock, ksblock_2)
226 7464 : IF (SIZE(ks_matrix, 1) > 1) THEN
227 : CALL dbcsr_get_block_p(matrix=ks_matrix(2, 1)%matrix, &
228 0 : row=iatom, col=jatom, block=ksblock_2, found=found)
229 : END IF
230 : CALL dbcsr_get_block_p(matrix=matrix_s(1)%matrix, &
231 7464 : row=iatom, col=jatom, block=sblock, found=found)
232 7464 : gmij = 0.5_dp*(gmcharge(iatom, 1) + gmcharge(jatom, 1))
233 129200 : ksblock = ksblock - gmij*sblock
234 7464 : IF (SIZE(ks_matrix, 1) > 1) ksblock_2 = ksblock_2 - gmij*sblock
235 9952 : IF (calculate_forces) THEN
236 24 : ikind = kind_of(iatom)
237 24 : atom_i = atom_of_kind(iatom)
238 24 : jkind = kind_of(jatom)
239 24 : atom_j = atom_of_kind(jatom)
240 24 : NULLIFY (pblock)
241 : CALL dbcsr_get_block_p(matrix=matrix_p(1)%matrix, &
242 24 : row=iatom, col=jatom, block=pblock, found=found)
243 96 : DO i = 1, 3
244 72 : NULLIFY (dsblock)
245 : CALL dbcsr_get_block_p(matrix=matrix_s(1 + i)%matrix, &
246 72 : row=iatom, col=jatom, block=dsblock, found=found)
247 696 : fi = -2.0_dp*gmij*SUM(pblock*dsblock)
248 72 : force(ikind)%rho_elec(i, atom_i) = force(ikind)%rho_elec(i, atom_i) + fi
249 72 : force(jkind)%rho_elec(i, atom_j) = force(jkind)%rho_elec(i, atom_j) - fi
250 96 : fij(i) = fi
251 : END DO
252 : END IF
253 : END DO
254 2488 : CALL dbcsr_iterator_stop(iter)
255 2488 : IF (calculate_forces .AND. SIZE(matrix_p) == 2) THEN
256 : CALL dbcsr_add(matrix_p(1)%matrix, matrix_p(2)%matrix, &
257 0 : alpha_scalar=1.0_dp, beta_scalar=-1.0_dp)
258 : END IF
259 : END IF
260 :
261 2488 : DEALLOCATE (gmcharge)
262 :
263 2488 : CALL timestop(handle)
264 :
265 7464 : END SUBROUTINE build_tb_coulomb_qmqm
266 :
267 : END MODULE qmmm_tb_coulomb
268 :
|
