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 Chemical shift calculation by dfpt
10 : !> Initialization of the nmr_env, creation of the special neighbor lists
11 : !> Perturbation Hamiltonians by application of the p and rxp oprtators to psi0
12 : !> Write output
13 : !> Deallocate everything
14 : !> \note
15 : !> The psi0 should be localized
16 : !> the Sebastiani method works within the assumption that the orbitals are
17 : !> completely contained in the simulation box
18 : !> \par History
19 : !> created 07-2005 [MI]
20 : !> \author MI
21 : ! **************************************************************************************************
22 : MODULE qs_linres_nmr_utils
23 :
24 : USE atomic_kind_types, ONLY: atomic_kind_type
25 : USE cell_types, ONLY: cell_type
26 : USE cp_control_types, ONLY: dft_control_type
27 : USE cp_log_handling, ONLY: cp_get_default_logger,&
28 : cp_logger_type
29 : USE cp_output_handling, ONLY: cp_p_file,&
30 : cp_print_key_finished_output,&
31 : cp_print_key_should_output,&
32 : cp_print_key_unit_nr
33 : USE cp_parser_methods, ONLY: parser_get_next_line,&
34 : parser_get_object
35 : USE cp_parser_types, ONLY: cp_parser_type,&
36 : parser_create,&
37 : parser_release
38 : USE cp_units, ONLY: cp_unit_from_cp2k,&
39 : cp_unit_to_cp2k
40 : USE input_section_types, ONLY: section_vals_get_subs_vals,&
41 : section_vals_type,&
42 : section_vals_val_get
43 : USE kinds, ONLY: default_path_length,&
44 : dp
45 : USE memory_utilities, ONLY: reallocate
46 : USE particle_types, ONLY: particle_type
47 : USE pw_env_types, ONLY: pw_env_type
48 : USE qs_environment_types, ONLY: get_qs_env,&
49 : qs_environment_type
50 : USE qs_linres_types, ONLY: linres_control_type,&
51 : nmr_env_type
52 : USE qs_matrix_pools, ONLY: qs_matrix_pools_type
53 : USE scf_control_types, ONLY: scf_control_type
54 : #include "./base/base_uses.f90"
55 :
56 : IMPLICIT NONE
57 :
58 : PRIVATE
59 : PUBLIC :: nmr_env_cleanup, nmr_env_init
60 :
61 : CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'qs_linres_nmr_utils'
62 :
63 : CONTAINS
64 :
65 : ! **************************************************************************************************
66 : !> \brief Initialize the nmr environment
67 : !> \param nmr_env ...
68 : !> \param qs_env ...
69 : !> \par History
70 : !> 07.2006 created [MI]
71 : !> \author MI
72 : ! **************************************************************************************************
73 160 : SUBROUTINE nmr_env_init(nmr_env, qs_env)
74 : !
75 : TYPE(nmr_env_type) :: nmr_env
76 : TYPE(qs_environment_type), POINTER :: qs_env
77 :
78 : CHARACTER(LEN=*), PARAMETER :: routineN = 'nmr_env_init'
79 :
80 : CHARACTER(LEN=default_path_length) :: nics_file_name
81 : INTEGER :: handle, ini, ir, j, n_mo(2), n_rep, nao, &
82 : nat_print, natom, nmoloc, nspins, &
83 : output_unit
84 160 : INTEGER, DIMENSION(:), POINTER :: bounds, list
85 : LOGICAL :: gapw
86 160 : TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
87 : TYPE(cell_type), POINTER :: cell
88 : TYPE(cp_logger_type), POINTER :: logger
89 : TYPE(dft_control_type), POINTER :: dft_control
90 : TYPE(linres_control_type), POINTER :: linres_control
91 160 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
92 : TYPE(pw_env_type), POINTER :: pw_env
93 : TYPE(qs_matrix_pools_type), POINTER :: mpools
94 : TYPE(scf_control_type), POINTER :: scf_control
95 : TYPE(section_vals_type), POINTER :: lr_section, nmr_section
96 :
97 : !
98 :
99 160 : CALL timeset(routineN, handle)
100 :
101 160 : NULLIFY (atomic_kind_set, cell, dft_control, linres_control, scf_control)
102 160 : NULLIFY (logger, mpools, nmr_section, particle_set)
103 160 : NULLIFY (pw_env)
104 :
105 : n_mo(1:2) = 0
106 160 : nao = 0
107 160 : nmoloc = 0
108 :
109 160 : logger => cp_get_default_logger()
110 160 : lr_section => section_vals_get_subs_vals(qs_env%input, "PROPERTIES%LINRES")
111 :
112 : output_unit = cp_print_key_unit_nr(logger, lr_section, "PRINT%PROGRAM_RUN_INFO", &
113 160 : extension=".linresLog")
114 :
115 160 : CALL nmr_env_cleanup(nmr_env)
116 :
117 160 : IF (output_unit > 0) THEN
118 80 : WRITE (output_unit, "(/,T20,A,/)") "*** Start NMR Chemical Shift Calculation ***"
119 80 : WRITE (output_unit, "(T10,A,/)") "Inizialization of the NMR environment"
120 : END IF
121 :
122 : ! If current_density or full_nmr different allocations are required
123 : nmr_section => section_vals_get_subs_vals(qs_env%input, &
124 160 : & "PROPERTIES%LINRES%NMR")
125 160 : CALL section_vals_val_get(nmr_section, "INTERPOLATE_SHIFT", l_val=nmr_env%interpolate_shift)
126 160 : CALL section_vals_val_get(nmr_section, "SHIFT_GAPW_RADIUS", r_val=nmr_env%shift_gapw_radius)
127 160 : CALL section_vals_val_get(nmr_section, "NICS", l_val=nmr_env%do_nics)
128 160 : IF (nmr_env%do_nics) THEN
129 : CALL section_vals_val_get(nmr_section, "NICS_FILE_NAME", &
130 6 : c_val=nics_file_name)
131 : BLOCK
132 : CHARACTER(LEN=2) :: label
133 : LOGICAL :: my_end
134 : TYPE(cp_parser_type) :: parser
135 6 : CALL parser_create(parser, nics_file_name)
136 6 : CALL parser_get_next_line(parser, 1)
137 6 : CALL parser_get_object(parser, nmr_env%n_nics)
138 18 : ALLOCATE (nmr_env%r_nics(3, nmr_env%n_nics))
139 6 : CALL parser_get_next_line(parser, 2)
140 24 : DO j = 1, nmr_env%n_nics
141 24 : CALL parser_get_object(parser, label)
142 24 : CALL parser_get_object(parser, nmr_env%r_nics(1, j))
143 24 : CALL parser_get_object(parser, nmr_env%r_nics(2, j))
144 24 : CALL parser_get_object(parser, nmr_env%r_nics(3, j))
145 24 : nmr_env%r_nics(1, j) = cp_unit_to_cp2k(nmr_env%r_nics(1, j), "angstrom")
146 24 : nmr_env%r_nics(2, j) = cp_unit_to_cp2k(nmr_env%r_nics(2, j), "angstrom")
147 24 : nmr_env%r_nics(3, j) = cp_unit_to_cp2k(nmr_env%r_nics(3, j), "angstrom")
148 24 : CALL parser_get_next_line(parser, 1, at_end=my_end)
149 24 : IF (my_end) EXIT
150 : END DO
151 24 : CALL parser_release(parser)
152 : END BLOCK
153 : END IF
154 :
155 : CALL get_qs_env(qs_env=qs_env, &
156 : atomic_kind_set=atomic_kind_set, &
157 : cell=cell, &
158 : dft_control=dft_control, &
159 : linres_control=linres_control, &
160 : mpools=mpools, &
161 : particle_set=particle_set, &
162 : pw_env=pw_env, &
163 160 : scf_control=scf_control)
164 : !
165 : ! Check if restat also psi0 should be restarted
166 : !IF(nmr_env%restart_nmr .AND. scf_control%density_guess/=restart_guess)THEN
167 : ! CPABORT("restart_nmr requires density_guess=restart")
168 : !ENDIF
169 : !
170 : ! check that the psi0 are localized and you have all the centers
171 160 : IF (.NOT. linres_control%localized_psi0) THEN
172 0 : CPWARN("To get NMR parameters within PBC you need localized zero order orbitals")
173 : END IF
174 160 : gapw = dft_control%qs_control%gapw
175 160 : nspins = dft_control%nspins
176 160 : natom = SIZE(particle_set, 1)
177 :
178 : !
179 : ! Conversion factors
180 : ! factor for the CHEMICAL SHIFTS: alpha^2 * ppm.
181 160 : nmr_env%shift_factor = (1.0_dp/137.03602_dp)**2*1.0E+6_dp/cell%deth
182 : ! factor for the CHEMICAL SHIFTS: alpha^2 * ppm.
183 160 : nmr_env%shift_factor_gapw = (1.0_dp/137.03602_dp)**2*1.0E+6_dp
184 : ! chi_factor = 1/4 * e^2/m * a_0 ^2
185 160 : nmr_env%chi_factor = 1.9727566E-29_dp/1.0E-30_dp ! -> displayed in 10^-30 J/T^2
186 : ! Factor to convert 10^-30 J/T^2 into ppm cgs = ppm cm^3/mol
187 : ! = 10^-30 * mu_0/4pi * N_A * 10^6 * 10^6 [one 10^6 for ppm, one for m^3 -> cm^3]
188 160 : nmr_env%chi_SI2ppmcgs = 6.022045_dp/1.0E+2_dp
189 : ! Chi to Shift: 10^-30 * 2/3 mu_0 / Omega * 1/ppm
190 : nmr_env%chi_SI2shiftppm = 1.0E-30_dp*8.37758041E-7_dp/ &
191 160 : (cp_unit_from_cp2k(cell%deth, "angstrom^3")*1.0E-30_dp)*1.0E+6_dp
192 :
193 160 : IF (output_unit > 0) THEN
194 80 : WRITE (output_unit, "(T2,A,T65,ES15.6)") "NMR| Shift gapw radius (a.u.) ", nmr_env%shift_gapw_radius
195 80 : IF (nmr_env%do_nics) THEN
196 3 : WRITE (output_unit, "(T2,A,T50,I5,A)") "NMR| NICS computed in ", nmr_env%n_nics, " additional points"
197 3 : WRITE (output_unit, "(T2,A,T60,A)") "NMR| NICS coordinates read on file ", TRIM(nics_file_name)
198 : END IF
199 80 : WRITE (output_unit, "(T2,A,T65,ES15.6)") "NMR| Shift factor (ppm)", nmr_env%shift_factor
200 80 : IF (gapw) THEN
201 43 : WRITE (output_unit, "(T2,A,T65,ES15.6)") "NMR| Shift factor gapw (ppm)", nmr_env%shift_factor_gapw
202 : END IF
203 80 : WRITE (output_unit, "(T2,A,T65,ES15.6)") "NMR| Chi factor (SI)", nmr_env%chi_factor
204 80 : WRITE (output_unit, "(T2,A,T65,ES15.6)") "NMR| Conversion Chi (ppm/cgs)", nmr_env%chi_SI2ppmcgs
205 80 : WRITE (output_unit, "(T2,A,T65,ES15.6)") "NMR| Conversion Chi to Shift", nmr_env%chi_SI2shiftppm
206 : END IF
207 :
208 480 : ALLOCATE (nmr_env%do_calc_cs_atom(natom))
209 654 : nmr_env%do_calc_cs_atom = 0
210 :
211 160 : IF (BTEST(cp_print_key_should_output(logger%iter_info, nmr_section,&
212 : & "PRINT%SHIELDING_TENSOR"), cp_p_file)) THEN
213 :
214 160 : NULLIFY (bounds, list)
215 160 : nat_print = 0
216 : CALL section_vals_val_get(nmr_section,&
217 : & "PRINT%SHIELDING_TENSOR%ATOMS_LU_BOUNDS", &
218 160 : i_vals=bounds)
219 160 : nat_print = bounds(2) - bounds(1) + 1
220 160 : IF (nat_print > 0) THEN
221 0 : ALLOCATE (nmr_env%cs_atom_list(nat_print))
222 0 : DO ir = 1, nat_print
223 0 : nmr_env%cs_atom_list(ir) = bounds(1) + (ir - 1)
224 0 : nmr_env%do_calc_cs_atom(bounds(1) + (ir - 1)) = 1
225 : END DO
226 : END IF
227 :
228 160 : IF (.NOT. ASSOCIATED(nmr_env%cs_atom_list)) THEN
229 : CALL section_vals_val_get(nmr_section, "PRINT%SHIELDING_TENSOR%ATOMS_LIST", &
230 160 : n_rep_val=n_rep)
231 160 : nat_print = 0
232 164 : DO ir = 1, n_rep
233 4 : NULLIFY (list)
234 : CALL section_vals_val_get(nmr_section, "PRINT%SHIELDING_TENSOR%ATOMS_LIST", &
235 4 : i_rep_val=ir, i_vals=list)
236 164 : IF (ASSOCIATED(list)) THEN
237 4 : CALL reallocate(nmr_env%cs_atom_list, 1, nat_print + SIZE(list))
238 10 : DO ini = 1, SIZE(list)
239 6 : nmr_env%cs_atom_list(ini + nat_print) = list(ini)
240 10 : nmr_env%do_calc_cs_atom(list(ini)) = 1
241 : END DO
242 4 : nat_print = nat_print + SIZE(list)
243 : END IF
244 : END DO ! ir
245 : END IF
246 :
247 160 : IF (.NOT. ASSOCIATED(nmr_env%cs_atom_list)) THEN
248 312 : ALLOCATE (nmr_env%cs_atom_list(natom))
249 638 : DO ir = 1, natom
250 638 : nmr_env%cs_atom_list(ir) = ir
251 : END DO
252 638 : nmr_env%do_calc_cs_atom = 1
253 : END IF
254 : !
255 : ! check the list
256 160 : CPASSERT(ASSOCIATED(nmr_env%cs_atom_list))
257 648 : DO ir = 1, SIZE(nmr_env%cs_atom_list, 1)
258 488 : IF (nmr_env%cs_atom_list(ir) .LT. 1 .OR. nmr_env%cs_atom_list(ir) .GT. natom) THEN
259 0 : CPABORT("Unknown atom(s)")
260 : END IF
261 2452 : DO j = 1, SIZE(nmr_env%cs_atom_list, 1)
262 1804 : IF (j .EQ. ir) CYCLE
263 1804 : IF (nmr_env%cs_atom_list(ir) .EQ. nmr_env%cs_atom_list(j)) THEN
264 0 : CPABORT("Duplicate atoms")
265 : END IF
266 : END DO
267 : END DO
268 : ELSE
269 0 : NULLIFY (nmr_env%cs_atom_list)
270 : END IF
271 :
272 160 : IF (output_unit > 0) THEN
273 80 : IF (ASSOCIATED(nmr_env%cs_atom_list)) THEN
274 80 : WRITE (output_unit, "(T2,A,T69,I5,A)") "NMR| Shielding tensor computed for ", &
275 160 : SIZE(nmr_env%cs_atom_list, 1), " atoms"
276 : ELSE
277 : WRITE (output_unit, "(T2,A,T50)")&
278 0 : & "NMR| Shielding tensor not computed at the atomic positions"
279 : END IF
280 : END IF
281 : !
282 : ! Initialize the chemical shift tensor
283 : ALLOCATE (nmr_env%chemical_shift(3, 3, natom), &
284 640 : nmr_env%chemical_shift_loc(3, 3, natom))
285 6582 : nmr_env%chemical_shift = 0.0_dp
286 6582 : nmr_env%chemical_shift_loc = 0.0_dp
287 160 : IF (nmr_env%do_nics) THEN
288 : ALLOCATE (nmr_env%chemical_shift_nics_loc(3, 3, nmr_env%n_nics), &
289 24 : nmr_env%chemical_shift_nics(3, 3, nmr_env%n_nics))
290 318 : nmr_env%chemical_shift_nics_loc = 0.0_dp
291 318 : nmr_env%chemical_shift_nics = 0.0_dp
292 : END IF
293 :
294 : CALL cp_print_key_finished_output(output_unit, logger, lr_section,&
295 160 : & "PRINT%PROGRAM_RUN_INFO")
296 :
297 160 : CALL timestop(handle)
298 :
299 160 : END SUBROUTINE nmr_env_init
300 :
301 : ! **************************************************************************************************
302 : !> \brief Deallocate the nmr environment
303 : !> \param nmr_env ...
304 : !> \par History
305 : !> 07.2005 created [MI]
306 : !> \author MI
307 : ! **************************************************************************************************
308 320 : SUBROUTINE nmr_env_cleanup(nmr_env)
309 :
310 : TYPE(nmr_env_type) :: nmr_env
311 :
312 320 : IF (ASSOCIATED(nmr_env%cs_atom_list)) THEN
313 160 : DEALLOCATE (nmr_env%cs_atom_list)
314 : END IF
315 320 : IF (ASSOCIATED(nmr_env%do_calc_cs_atom)) THEN
316 160 : DEALLOCATE (nmr_env%do_calc_cs_atom)
317 : END IF
318 : !chemical_shift
319 320 : IF (ASSOCIATED(nmr_env%chemical_shift)) THEN
320 160 : DEALLOCATE (nmr_env%chemical_shift)
321 : END IF
322 320 : IF (ASSOCIATED(nmr_env%chemical_shift_loc)) THEN
323 160 : DEALLOCATE (nmr_env%chemical_shift_loc)
324 : END IF
325 : ! nics
326 320 : IF (ASSOCIATED(nmr_env%r_nics)) THEN
327 6 : DEALLOCATE (nmr_env%r_nics)
328 : END IF
329 320 : IF (ASSOCIATED(nmr_env%chemical_shift_nics)) THEN
330 6 : DEALLOCATE (nmr_env%chemical_shift_nics)
331 : END IF
332 320 : IF (ASSOCIATED(nmr_env%chemical_shift_nics_loc)) THEN
333 6 : DEALLOCATE (nmr_env%chemical_shift_nics_loc)
334 : END IF
335 :
336 320 : END SUBROUTINE nmr_env_cleanup
337 :
338 : END MODULE qs_linres_nmr_utils
|
