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 : MODULE tip_scan_methods
9 : USE cp_control_types, ONLY: dft_control_type
10 : USE cp_files, ONLY: close_file,&
11 : open_file
12 : USE cp_log_handling, ONLY: cp_get_default_logger,&
13 : cp_logger_get_default_io_unit,&
14 : cp_logger_type
15 : USE cp_output_handling, ONLY: silent_print_level
16 : USE cp_realspace_grid_cube, ONLY: cp_cube_to_pw
17 : USE input_section_types, ONLY: section_vals_type,&
18 : section_vals_val_get
19 : USE kinds, ONLY: default_string_length,&
20 : dp
21 : USE message_passing, ONLY: mp_para_env_type
22 : USE pw_env_types, ONLY: pw_env_get,&
23 : pw_env_type
24 : USE pw_grid_types, ONLY: pw_grid_type
25 : USE pw_grids, ONLY: pw_grid_create,&
26 : pw_grid_release
27 : USE pw_methods, ONLY: pw_axpy,&
28 : pw_multiply_with,&
29 : pw_structure_factor,&
30 : pw_transfer,&
31 : pw_zero
32 : USE pw_pool_types, ONLY: pw_pool_type
33 : USE pw_types, ONLY: pw_c1d_gs_type,&
34 : pw_r3d_rs_type
35 : USE qs_environment_types, ONLY: get_qs_env,&
36 : qs_environment_type
37 : USE qs_ks_types, ONLY: qs_ks_env_type,&
38 : set_ks_env
39 : USE qs_mo_types, ONLY: deallocate_mo_set,&
40 : duplicate_mo_set,&
41 : mo_set_type,&
42 : reassign_allocated_mos
43 : USE qs_scf, ONLY: scf
44 : USE tip_scan_types, ONLY: read_scanning_section,&
45 : release_scanning_type,&
46 : scanning_type
47 : #include "./base/base_uses.f90"
48 :
49 : IMPLICIT NONE
50 :
51 : PRIVATE
52 :
53 : CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'tip_scan_methods'
54 :
55 : PUBLIC :: tip_scanning
56 :
57 : ! **************************************************************************************************
58 :
59 : CONTAINS
60 :
61 : ! **************************************************************************************************
62 : !> \brief Perform tip scanning calculation.
63 : !> \param qs_env Quickstep environment
64 : !> input_section Tip Scan Section
65 : !> \param input_section ...
66 : !> \par History
67 : !> * 05.2021 created [JGH]
68 : ! **************************************************************************************************
69 0 : SUBROUTINE tip_scanning(qs_env, input_section)
70 : TYPE(qs_environment_type), POINTER :: qs_env
71 : TYPE(section_vals_type), POINTER :: input_section
72 :
73 : CHARACTER(LEN=*), PARAMETER :: routineN = 'tip_scanning'
74 :
75 : CHARACTER(LEN=default_string_length) :: cname
76 : INTEGER :: handle, iounit, iscan, iset, nscan, &
77 : nset, plevel, tsteps
78 : LOGICAL :: do_tip_scan, expot, scf_converged
79 : REAL(KIND=dp), DIMENSION(3) :: rpos
80 : TYPE(cp_logger_type), POINTER :: logger
81 : TYPE(dft_control_type), POINTER :: dft_control
82 0 : TYPE(mo_set_type), ALLOCATABLE, DIMENSION(:) :: mos_ref
83 0 : TYPE(mo_set_type), DIMENSION(:), POINTER :: mos
84 : TYPE(pw_c1d_gs_type) :: sf, vref
85 : TYPE(pw_env_type), POINTER :: pw_env
86 : TYPE(pw_pool_type), POINTER :: auxbas_pw_pool
87 : TYPE(pw_r3d_rs_type), POINTER :: vee, vtip
88 : TYPE(qs_ks_env_type), POINTER :: ks_env
89 0 : TYPE(scanning_type) :: scan_env
90 :
91 0 : CALL timeset(routineN, handle)
92 :
93 0 : NULLIFY (logger)
94 0 : logger => cp_get_default_logger()
95 :
96 0 : CALL section_vals_val_get(input_section, "_SECTION_PARAMETERS_", l_val=do_tip_scan)
97 0 : IF (do_tip_scan) THEN
98 0 : iounit = cp_logger_get_default_io_unit(logger)
99 0 : cname = logger%iter_info%project_name
100 0 : logger%iter_info%project_name = TRIM(logger%iter_info%project_name)//"+TIP_SCAN"
101 0 : plevel = logger%iter_info%print_level
102 0 : logger%iter_info%print_level = silent_print_level
103 :
104 0 : IF (iounit > 0) THEN
105 0 : WRITE (iounit, "(T2,A)") "TIP SCAN| Perform a Tip Scanning Calculation"
106 : END IF
107 :
108 : ! read the input section
109 0 : CALL read_scanning_section(scan_env, input_section)
110 : ! read tip potential file
111 0 : CALL read_tip_file(qs_env, scan_env)
112 :
113 : CALL get_qs_env(qs_env, ks_env=ks_env, pw_env=pw_env, &
114 0 : dft_control=dft_control)
115 0 : expot = dft_control%apply_external_potential
116 0 : dft_control%apply_external_potential = .TRUE.
117 0 : IF (expot) THEN
118 : ! save external potential
119 0 : CALL get_qs_env(qs_env, vee=vee)
120 : END IF
121 :
122 : ! scratch memory for tip potentials and structure factor
123 0 : CALL pw_env_get(pw_env, auxbas_pw_pool=auxbas_pw_pool)
124 : NULLIFY (vtip)
125 0 : ALLOCATE (vtip)
126 0 : CALL auxbas_pw_pool%create_pw(vtip)
127 0 : CALL pw_zero(vtip)
128 0 : CALL auxbas_pw_pool%create_pw(vref)
129 0 : CALL pw_zero(vref)
130 0 : CALL auxbas_pw_pool%create_pw(sf)
131 :
132 : ! get the reference tip potential and store it in reciprocal space (vref)
133 0 : CALL pw_transfer(scan_env%tip_pw_g, vref)
134 :
135 : ! store reference MOs
136 0 : CALL get_qs_env(qs_env, mos=mos)
137 0 : nset = SIZE(mos)
138 0 : ALLOCATE (mos_ref(nset))
139 0 : DO iset = 1, nset
140 0 : CALL duplicate_mo_set(mos_ref(iset), mos(iset))
141 : END DO
142 :
143 0 : nscan = scan_env%num_scan_points
144 0 : IF (iounit > 0) THEN
145 0 : WRITE (iounit, "(T2,A,T74,I7)") "TIP SCAN| Number of scanning points ", nscan
146 0 : WRITE (iounit, "(T2,A)") "TIP SCAN| Start scanning ..."
147 : END IF
148 :
149 0 : DO iscan = 1, nscan
150 0 : IF (iounit > 0) THEN
151 0 : WRITE (iounit, "(T2,A,I7)", advance="NO") "TIP SCAN| Scan point ", iscan
152 : END IF
153 :
154 : ! shift the reference tip potential
155 0 : rpos(1:3) = scan_env%tip_pos(1:3, iscan) - scan_env%ref_point(1:3)
156 0 : CALL shift_tip_potential(vref, sf, vtip, rpos)
157 : ! set the external potential
158 0 : IF (ASSOCIATED(vee)) THEN
159 0 : CALL pw_axpy(vee, vtip, alpha=1.0_dp)
160 : END IF
161 0 : CALL set_ks_env(ks_env, vee=vtip)
162 :
163 : ! reset MOs
164 0 : CALL get_qs_env(qs_env, mos=mos)
165 0 : DO iset = 1, nset
166 0 : CALL reassign_allocated_mos(mos(iset), mos_ref(iset))
167 : END DO
168 :
169 : ! Calculate electronic structure
170 0 : CALL scf(qs_env, has_converged=scf_converged, total_scf_steps=tsteps)
171 :
172 0 : IF (iounit > 0) THEN
173 0 : IF (scf_converged) THEN
174 0 : WRITE (iounit, "(T25,A,I4,A)") "SCF converged in ", tsteps, " steps"
175 : ELSE
176 0 : WRITE (iounit, "(T31,A)") "SCF did not converge!"
177 : END IF
178 : END IF
179 : END DO
180 0 : CALL release_scanning_type(scan_env)
181 :
182 0 : IF (iounit > 0) THEN
183 0 : WRITE (iounit, "(T2,A)") "TIP SCAN| ... end scanning"
184 : END IF
185 0 : dft_control%apply_external_potential = expot
186 0 : IF (expot) THEN
187 : ! restore vee
188 0 : CALL set_ks_env(ks_env, vee=vee)
189 : ELSE
190 0 : NULLIFY (vee)
191 0 : CALL set_ks_env(ks_env, vee=vee)
192 : END IF
193 0 : CALL auxbas_pw_pool%give_back_pw(vtip)
194 0 : CALL auxbas_pw_pool%give_back_pw(vref)
195 0 : CALL auxbas_pw_pool%give_back_pw(sf)
196 0 : DEALLOCATE (vtip)
197 :
198 0 : logger%iter_info%print_level = plevel
199 0 : logger%iter_info%project_name = cname
200 :
201 : ! reset MOs
202 0 : CALL get_qs_env(qs_env, mos=mos)
203 0 : DO iset = 1, nset
204 0 : CALL reassign_allocated_mos(mos(iset), mos_ref(iset))
205 0 : CALL deallocate_mo_set(mos_ref(iset))
206 : END DO
207 0 : DEALLOCATE (mos_ref)
208 : END IF
209 :
210 0 : CALL timestop(handle)
211 :
212 0 : END SUBROUTINE tip_scanning
213 :
214 : ! **************************************************************************************************
215 : !> \brief Shift tip potential in reciprocal space
216 : !> \param vref ...
217 : !> \param sf ...
218 : !> \param vtip ...
219 : !> \param rpos ...
220 : !> \par History
221 : !> * 05.2021 created [JGH]
222 : ! **************************************************************************************************
223 0 : SUBROUTINE shift_tip_potential(vref, sf, vtip, rpos)
224 :
225 : TYPE(pw_c1d_gs_type), INTENT(INOUT) :: vref, sf
226 : TYPE(pw_r3d_rs_type), INTENT(INOUT) :: vtip
227 : REAL(KIND=dp), DIMENSION(3), INTENT(IN) :: rpos
228 :
229 : CHARACTER(LEN=*), PARAMETER :: routineN = 'shift_tip_potential'
230 :
231 : INTEGER :: handle
232 :
233 0 : CALL timeset(routineN, handle)
234 :
235 0 : CALL pw_structure_factor(sf, rpos)
236 0 : CALL pw_multiply_with(sf, vref)
237 0 : CALL pw_transfer(sf, vtip)
238 :
239 0 : CALL timestop(handle)
240 :
241 0 : END SUBROUTINE shift_tip_potential
242 :
243 : ! **************************************************************************************************
244 : !> \brief Read tip potential from cube file. Allow any spacing and cell size
245 : !> \param qs_env ...
246 : !> \param scan_env ...
247 : !> \par History
248 : !> * 05.2021 created [JGH]
249 : ! **************************************************************************************************
250 0 : SUBROUTINE read_tip_file(qs_env, scan_env)
251 : TYPE(qs_environment_type), POINTER :: qs_env
252 : TYPE(scanning_type), INTENT(INOUT) :: scan_env
253 :
254 : CHARACTER(LEN=*), PARAMETER :: routineN = 'read_tip_file'
255 :
256 : INTEGER :: extunit, handle, i, nat
257 : INTEGER, DIMENSION(3) :: npts
258 : REAL(KIND=dp) :: scaling
259 : REAL(KIND=dp), DIMENSION(3) :: rdum
260 : REAL(KIND=dp), DIMENSION(3, 3) :: dcell
261 : TYPE(mp_para_env_type), POINTER :: para_env
262 : TYPE(pw_grid_type), POINTER :: pw_grid
263 :
264 0 : CALL timeset(routineN, handle)
265 :
266 0 : CALL get_qs_env(qs_env, para_env=para_env)
267 :
268 0 : IF (para_env%is_source()) THEN
269 : CALL open_file(file_name=scan_env%tip_cube_file, &
270 : file_status="OLD", &
271 : file_form="FORMATTED", &
272 : file_action="READ", &
273 0 : unit_number=extunit)
274 : !skip header comments
275 0 : DO i = 1, 2
276 0 : READ (extunit, *)
277 : END DO
278 0 : READ (extunit, *) nat, rdum
279 0 : DO i = 1, 3
280 0 : READ (extunit, *) npts(i), dcell(i, 1:3)
281 0 : dcell(i, 1:3) = npts(i)*dcell(i, 1:3)
282 : END DO
283 0 : CALL close_file(unit_number=extunit)
284 : END IF
285 :
286 0 : CALL para_env%bcast(npts)
287 0 : CALL para_env%bcast(dcell)
288 :
289 0 : NULLIFY (pw_grid)
290 0 : CALL pw_grid_create(pw_grid, para_env, dcell, npts=npts)
291 0 : CALL scan_env%tip_pw_r%create(pw_grid)
292 : !deb
293 0 : scaling = 0.1_dp
294 : !deb
295 0 : CALL cp_cube_to_pw(scan_env%tip_pw_r, scan_env%tip_cube_file, scaling, silent=.TRUE.)
296 0 : CALL scan_env%tip_pw_g%create(pw_grid)
297 0 : CALL pw_transfer(scan_env%tip_pw_r, scan_env%tip_pw_g)
298 0 : CALL pw_grid_release(pw_grid)
299 :
300 0 : CALL timestop(handle)
301 :
302 0 : END SUBROUTINE read_tip_file
303 :
304 : END MODULE tip_scan_methods
|
