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