Line data Source code
1 : !--------------------------------------------------------------------------------------------------!
2 : ! CP2K: A general program to perform molecular dynamics simulations !
3 : ! Copyright 2000-2026 CP2K developers group <https://cp2k.org> !
4 : ! !
5 : ! SPDX-License-Identifier: GPL-2.0-or-later !
6 : !--------------------------------------------------------------------------------------------------!
7 : MODULE fde_main
8 : USE bibliography, ONLY: Schreder2024_3,&
9 : cite_reference
10 : USE cp_control_types, ONLY: dft_control_type,&
11 : fde_control_type
12 : USE cp_files, ONLY: join_paths
13 : USE cp_log_handling, ONLY: cp_get_default_logger,&
14 : cp_logger_get_default_io_unit,&
15 : cp_logger_type
16 : USE fde_methods, ONLY: fde_build_v_emb,&
17 : fde_read_pw,&
18 : fde_read_results,&
19 : fde_spawn,&
20 : fde_v_es,&
21 : fde_write_pw
22 : USE kinds, ONLY: default_path_length,&
23 : dp
24 : USE message_passing, ONLY: mp_para_env_type
25 : USE pw_env_types, ONLY: pw_env_get,&
26 : pw_env_type
27 : USE pw_methods, ONLY: pw_axpy,&
28 : pw_integral_ab,&
29 : pw_integrate_function,&
30 : pw_scale,&
31 : pw_zero
32 : USE pw_pool_types, ONLY: pw_pool_type
33 : USE pw_types, ONLY: pw_r3d_rs_type
34 : USE qs_energy_types, ONLY: qs_energy_type
35 : USE qs_environment_types, ONLY: get_qs_env,&
36 : qs_environment_type
37 : USE qs_external_potential, ONLY: interpolate_external_potential
38 : #include "./base/base_uses.f90"
39 :
40 : IMPLICIT NONE
41 : PRIVATE
42 :
43 : CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'fde_main'
44 :
45 : PUBLIC :: fde_run
46 :
47 : CONTAINS
48 :
49 : ! **************************************************************************************************
50 : !> \brief Run the FDE self-consistency loop
51 : !> \param qs_env ..
52 : ! **************************************************************************************************
53 2 : SUBROUTINE fde_run(qs_env)
54 : TYPE(qs_environment_type), POINTER :: qs_env
55 :
56 : CHARACTER(len=*), PARAMETER :: routineN = 'fde_run'
57 :
58 : CHARACTER(len=default_path_length) :: grid_path, input, res_path, rho_path, &
59 : v_path
60 : INTEGER :: handle, i_a, iounit, iter, max_iter, n_a
61 : LOGICAL :: converged
62 : REAL(KIND=dp) :: drho, e_emb, e_env, e_es_e, e_es_int, &
63 : e_es_n, e_kin_na, e_n, e_sub, &
64 : e_sub_bare, e_tot, e_xc_na, n_el, &
65 : n_el_pw, phi_at_atom
66 2 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: z
67 2 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: r
68 : TYPE(cp_logger_type), POINTER :: logger
69 : TYPE(dft_control_type), POINTER :: dft_control
70 : TYPE(fde_control_type), POINTER :: fde_control
71 : TYPE(mp_para_env_type), POINTER :: para_env
72 : TYPE(pw_env_type), POINTER :: pw_env
73 : TYPE(pw_pool_type), POINTER :: auxbas_pw_pool
74 : TYPE(pw_r3d_rs_type) :: rho, rho_new, v_emb
75 : TYPE(pw_r3d_rs_type), POINTER :: phi_es_ptr
76 : TYPE(pw_r3d_rs_type), TARGET :: phi_es
77 : TYPE(qs_energy_type), POINTER :: qs_energy
78 :
79 2 : CALL timeset(routineN, handle)
80 :
81 2 : CALL cite_reference(Schreder2024_3)
82 :
83 2 : logger => cp_get_default_logger()
84 2 : iounit = cp_logger_get_default_io_unit(logger)
85 :
86 2 : NULLIFY (dft_control, fde_control, para_env, pw_env, auxbas_pw_pool, qs_energy)
87 : CALL get_qs_env(qs_env, dft_control=dft_control, para_env=para_env, pw_env=pw_env, &
88 2 : energy=qs_energy)
89 2 : fde_control => dft_control%fde_control
90 2 : CPASSERT(ASSOCIATED(fde_control))
91 :
92 2 : CALL pw_env_get(pw_env, auxbas_pw_pool=auxbas_pw_pool)
93 :
94 : ! socket files
95 2 : v_path = join_paths(fde_control%workdir, fde_control%pot_file)
96 2 : grid_path = join_paths(fde_control%workdir, fde_control%grid_file)
97 2 : rho_path = join_paths(fde_control%workdir, fde_control%rho_file)
98 2 : res_path = join_paths(fde_control%workdir, fde_control%result_file)
99 2 : input = fde_control%molcas_input_template
100 :
101 : ! bare environment DFT energy
102 2 : e_env = qs_energy%total
103 :
104 : ! initial guess = zero
105 2 : CALL auxbas_pw_pool%create_pw(rho)
106 2 : CALL pw_zero(rho)
107 2 : CALL auxbas_pw_pool%create_pw(rho_new)
108 :
109 2 : max_iter = 1
110 2 : IF (fde_control%freeze_and_thaw) max_iter = fde_control%max_ft_iter
111 :
112 2 : IF (iounit > 0) THEN
113 1 : WRITE (iounit, '(/,T2,A)') REPEAT("=", 78)
114 1 : WRITE (iounit, '(T2,A)') "FDE: Frozen-density embedding"
115 1 : WRITE (iounit, '(T2,A)') REPEAT("=", 78)
116 : END IF
117 :
118 2 : converged = .FALSE.
119 4 : DO iter = 1, max_iter
120 4 : CALL fde_build_v_emb(qs_env, rho, v_emb, e_xc_na, e_kin_na)
121 4 : CALL fde_write_pw(v_emb, TRIM(v_path), write_value=.TRUE.)
122 4 : CALL fde_write_pw(v_emb, TRIM(grid_path), write_value=.FALSE.)
123 :
124 4 : CALL fde_spawn(fde_control, TRIM(input), para_env)
125 :
126 4 : CALL fde_read_pw(rho_new, TRIM(rho_path))
127 4 : CALL fde_read_results(TRIM(res_path), para_env, e_sub, e_emb, e_n, n_el, n_a, z, r)
128 :
129 4 : n_el_pw = pw_integrate_function(rho_new)
130 4 : IF (n_el_pw > 1.0e-10_dp .AND. n_el > 0.0_dp) THEN
131 4 : CALL pw_scale(rho_new, n_el/n_el_pw)
132 : END IF
133 :
134 4 : drho = fde_drho(rho, rho_new, para_env)
135 :
136 4 : CALL fde_v_es(qs_env, phi_es)
137 4 : phi_es_ptr => phi_es
138 4 : e_es_e = pw_integral_ab(rho_new, phi_es)
139 :
140 4 : e_es_n = 0.0_dp
141 8 : DO i_a = 1, n_a
142 4 : CALL interpolate_external_potential(r(:, i_a), phi_es_ptr, func=phi_at_atom)
143 8 : e_es_n = e_es_n - z(i_a)*phi_at_atom
144 : END DO
145 4 : e_es_int = e_es_e + e_es_n
146 4 : CALL auxbas_pw_pool%give_back_pw(phi_es)
147 4 : DEALLOCATE (z, r)
148 :
149 4 : e_sub_bare = e_sub - e_emb
150 4 : e_tot = e_sub_bare + e_env + e_es_int + e_xc_na + e_kin_na
151 :
152 4 : IF (iounit > 0) THEN
153 2 : WRITE (iounit, '(/,T2,A,I0)') "FDE freeze-and-thaw iteration ", iter
154 2 : WRITE (iounit, '(T4,A,F20.10)') "E(sub, embedded) = ", e_sub
155 2 : WRITE (iounit, '(T4,A,F20.10)') "E(emb) = ", e_emb
156 2 : WRITE (iounit, '(T4,A,F20.10)') "E(sub) = ", e_sub_bare
157 2 : WRITE (iounit, '(T4,A,F20.10)') "E(env) = ", e_env
158 2 : WRITE (iounit, '(T4,A,F20.10)') "E(es) = ", e_es_int
159 2 : WRITE (iounit, '(T6,A,F20.10)') "E(elec) = ", e_es_e
160 2 : WRITE (iounit, '(T6,A,F20.10)') "E(nuc) = ", e_es_n
161 2 : WRITE (iounit, '(T4,A,F20.10)') "E(XC, nonadditive) = ", e_xc_na
162 2 : WRITE (iounit, '(T4,A,F20.10)') "E(kin, nonadditive) = ", e_kin_na
163 2 : WRITE (iounit, '(T4,A,F20.10)') "E(FDE) = ", e_tot
164 2 : WRITE (iounit, '(T4,A,F14.6,A,F14.6)') "N(sub, e): grid = ", n_el_pw, &
165 4 : " reported = ", n_el
166 2 : WRITE (iounit, '(T4,A,ES14.6)') "int|rho(sub, new) - rho(sub, old)|dr = ", drho
167 : END IF
168 :
169 4 : CALL auxbas_pw_pool%give_back_pw(v_emb)
170 :
171 4 : IF (fde_control%freeze_and_thaw) THEN
172 4 : IF (drho < fde_control%eps_ft_density) THEN
173 : converged = .TRUE.
174 : ELSE
175 2 : CALL pw_scale(rho, 1.0_dp - fde_control%density_mixing)
176 2 : CALL pw_axpy(rho_new, rho, fde_control%density_mixing)
177 : END IF
178 : ELSE
179 : converged = .TRUE.
180 : END IF
181 :
182 8 : IF (converged) EXIT
183 :
184 : END DO
185 :
186 2 : IF (iounit > 0) THEN
187 1 : WRITE (iounit, '(/,T2,A)') REPEAT("=", 78)
188 1 : IF (fde_control%freeze_and_thaw .AND. .NOT. converged) THEN
189 0 : WRITE (iounit, '(T2,A)') "FDE WARNING: freeze-and-thaw did not converge"
190 : END IF
191 1 : WRITE (iounit, '(T2,A,F20.10)') "FDE final total energy = ", e_tot
192 1 : WRITE (iounit, '(T2,A)') REPEAT("=", 78)
193 : END IF
194 :
195 2 : qs_energy%total = e_tot
196 :
197 2 : CALL auxbas_pw_pool%give_back_pw(rho)
198 2 : CALL auxbas_pw_pool%give_back_pw(rho_new)
199 :
200 2 : CALL timestop(handle)
201 :
202 4 : END SUBROUTINE fde_run
203 :
204 : ! **************************************************************************************************
205 : !> \brief Integrated density difference: int |rho_new - rho_old| dr
206 : !> \param rho_a ...
207 : !> \param rho_b ...
208 : !> \param para_env ...
209 : !> \return ...
210 : ! **************************************************************************************************
211 4 : FUNCTION fde_drho(rho_a, rho_b, para_env) RESULT(drho)
212 : TYPE(pw_r3d_rs_type), INTENT(IN) :: rho_a, rho_b
213 : TYPE(mp_para_env_type), INTENT(IN) :: para_env
214 : REAL(KIND=dp) :: drho
215 :
216 13036 : drho = SUM(ABS(rho_b%array - rho_a%array))*rho_b%pw_grid%dvol
217 4 : CALL para_env%sum(drho)
218 :
219 4 : END FUNCTION fde_drho
220 :
221 : END MODULE fde_main
|
