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 : !> \brief Routines to handle an external density
9 : !> The external density can be generic and is provided by user input
10 : !> \author D. Varsano
11 : ! **************************************************************************************************
12 : MODULE qs_external_density
13 : USE cell_types, ONLY: cell_type
14 : USE cp_control_types, ONLY: dft_control_type
15 : USE cp_files, ONLY: close_file,&
16 : open_file
17 : USE gaussian_gridlevels, ONLY: gridlevel_info_type
18 : USE input_section_types, ONLY: section_vals_get_subs_vals,&
19 : section_vals_type,&
20 : section_vals_val_get
21 : USE kinds, ONLY: default_string_length,&
22 : dp
23 : USE pw_env_types, ONLY: pw_env_get,&
24 : pw_env_type
25 : USE pw_methods, ONLY: pw_integrate_function
26 : USE pw_types, ONLY: pw_c1d_gs_type,&
27 : pw_r3d_rs_type
28 : USE qs_environment_types, ONLY: get_qs_env,&
29 : qs_environment_type
30 : USE qs_rho_types, ONLY: qs_rho_get,&
31 : qs_rho_type
32 : USE realspace_grid_types, ONLY: realspace_grid_desc_p_type,&
33 : realspace_grid_type,&
34 : rs_grid_create,&
35 : rs_grid_release,&
36 : rs_grid_zero
37 : USE rs_pw_interface, ONLY: density_rs2pw
38 : #include "./base/base_uses.f90"
39 :
40 : IMPLICIT NONE
41 :
42 : PRIVATE
43 :
44 : CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'qs_external_density'
45 :
46 : PUBLIC :: external_read_density
47 :
48 : CONTAINS
49 :
50 : ! **************************************************************************************************
51 : !> \brief Computes the external density on the grid
52 : !> \param qs_env ...
53 : !> \date 03.2011
54 : !> \author D. Varsano
55 : ! **************************************************************************************************
56 10250 : SUBROUTINE external_read_density(qs_env)
57 :
58 : TYPE(qs_environment_type), POINTER :: qs_env
59 :
60 : CHARACTER(len=*), PARAMETER :: routineN = 'external_read_density'
61 :
62 : CHARACTER(LEN=default_string_length) :: filename
63 : INTEGER :: extunit, handle, i, igrid_level, j, k, &
64 : nat, ndum, tag
65 : INTEGER, DIMENSION(3) :: lbounds, lbounds_local, npoints, &
66 : npoints_local, ubounds, ubounds_local
67 10250 : REAL(kind=dp), ALLOCATABLE, DIMENSION(:) :: buffer
68 : REAL(kind=dp), DIMENSION(3) :: dr, rdum
69 10250 : REAL(KIND=dp), DIMENSION(:), POINTER :: tot_rho_r_ext
70 : TYPE(cell_type), POINTER :: cell
71 : TYPE(dft_control_type), POINTER :: dft_control
72 : TYPE(gridlevel_info_type), POINTER :: gridlevel_info
73 10250 : TYPE(pw_c1d_gs_type), DIMENSION(:), POINTER :: rho_ext_g
74 : TYPE(pw_env_type), POINTER :: pw_env
75 10250 : TYPE(pw_r3d_rs_type), DIMENSION(:), POINTER :: rho_ext_r
76 : TYPE(qs_rho_type), POINTER :: rho_external
77 : TYPE(realspace_grid_desc_p_type), DIMENSION(:), &
78 10250 : POINTER :: rs_descs
79 : TYPE(realspace_grid_type), ALLOCATABLE, &
80 10250 : DIMENSION(:) :: rs_rho_ext
81 : TYPE(section_vals_type), POINTER :: ext_den_section, input
82 :
83 10250 : CALL timeset(routineN, handle)
84 10250 : NULLIFY (cell, input, ext_den_section, rs_descs, dft_control)
85 10250 : NULLIFY (rho_ext_r, rho_ext_g, tot_rho_r_ext)
86 :
87 : CALL get_qs_env(qs_env, &
88 : cell=cell, &
89 : rho_external=rho_external, &
90 : input=input, &
91 : pw_env=pw_env, &
92 10250 : dft_control=dft_control)
93 :
94 10250 : IF (dft_control%apply_external_density) THEN
95 : CALL qs_rho_get(rho_external, &
96 : rho_r=rho_ext_r, &
97 : rho_g=rho_ext_g, &
98 0 : tot_rho_r=tot_rho_r_ext)
99 :
100 0 : gridlevel_info => pw_env%gridlevel_info
101 :
102 0 : CALL pw_env_get(pw_env, rs_descs=rs_descs)
103 :
104 0 : ALLOCATE (rs_rho_ext(gridlevel_info%ngrid_levels))
105 :
106 0 : DO igrid_level = 1, gridlevel_info%ngrid_levels
107 : CALL rs_grid_create(rs_rho_ext(igrid_level), &
108 0 : rs_descs(igrid_level)%rs_desc)
109 0 : CALL rs_grid_zero(rs_rho_ext(igrid_level))
110 : END DO
111 :
112 0 : igrid_level = igrid_level - 1
113 :
114 0 : ext_den_section => section_vals_get_subs_vals(input, "DFT%EXTERNAL_DENSITY")
115 0 : CALL section_vals_val_get(ext_den_section, "FILE_DENSITY", c_val=filename)
116 :
117 0 : tag = 1
118 : ASSOCIATE (gid => rho_ext_r(1)%pw_grid%para%group, my_rank => rho_ext_r(1)%pw_grid%para%my_pos, &
119 : num_pe => rho_ext_r(1)%pw_grid%para%group_size)
120 :
121 0 : IF (dft_control%read_external_density) THEN
122 :
123 0 : DO i = 1, 3
124 0 : dr(i) = rs_descs(igrid_level)%rs_desc%dh(i, i)
125 : END DO
126 0 : npoints = rs_descs(igrid_level)%rs_desc%npts
127 0 : lbounds = rs_descs(igrid_level)%rs_desc%lb
128 0 : ubounds = rs_descs(igrid_level)%rs_desc%ub
129 :
130 0 : npoints_local = rho_ext_r(1)%pw_grid%npts_local
131 0 : lbounds_local = rho_ext_r(1)%pw_grid%bounds_local(1, :)
132 0 : ubounds_local = rho_ext_r(1)%pw_grid%bounds_local(2, :)
133 :
134 0 : ALLOCATE (buffer(lbounds_local(3):ubounds_local(3)))
135 :
136 0 : IF (my_rank == 0) THEN
137 0 : WRITE (*, FMT="(/,/,T2,A)") "INITIALIZING ZMP CONSTRAINED DENSITY METHOD"
138 0 : WRITE (*, FMT="(/,(T3,A,T51,A30))") "ZMP| Reading the target density: ", filename
139 :
140 : CALL open_file(file_name=filename, &
141 : file_status="OLD", &
142 : file_form="FORMATTED", &
143 : file_action="READ", &
144 0 : unit_number=extunit)
145 :
146 0 : DO i = 1, 2
147 0 : READ (extunit, *)
148 : END DO
149 0 : READ (extunit, *) nat, rdum
150 0 : DO i = 1, 3
151 0 : READ (extunit, *) ndum, rdum
152 0 : IF (ndum /= npoints(i) .OR. (ABS(rdum(i) - dr(i)) > 1e-4)) THEN
153 0 : WRITE (*, *) "ZMP | ERROR! | CUBE FILE NOT COINCIDENT WITH INTERNAL GRID ", i
154 0 : WRITE (*, *) "ZMP | ", ndum, " DIFFERS FROM ", npoints(i)
155 0 : WRITE (*, *) "ZMP | ", rdum, " DIFFERS FROM ", dr(i)
156 : END IF
157 : END DO
158 0 : DO i = 1, nat
159 0 : READ (extunit, *)
160 : END DO
161 : END IF
162 :
163 0 : DO i = lbounds(1), ubounds(1)
164 0 : DO j = lbounds(2), ubounds(2)
165 0 : IF (my_rank .EQ. 0) THEN
166 0 : READ (extunit, *) (buffer(k), k=lbounds(3), ubounds(3))
167 : END IF
168 0 : CALL gid%bcast(buffer(lbounds(3):ubounds(3)), 0)
169 :
170 : IF ((lbounds_local(1) .LE. i) .AND. (i .LE. ubounds_local(1)) .AND. (lbounds_local(2) .LE. j) &
171 0 : .AND. (j .LE. ubounds_local(2))) THEN
172 0 : rs_rho_ext(igrid_level)%r(i, j, lbounds(3):ubounds(3)) = buffer(lbounds(3):ubounds(3))
173 : END IF
174 :
175 : END DO
176 : END DO
177 0 : IF (my_rank == 0) CALL close_file(unit_number=extunit)
178 : END IF
179 :
180 0 : CALL density_rs2pw(pw_env, rs_rho_ext, rho=rho_ext_r(1), rho_gspace=rho_ext_g(1))
181 0 : DO igrid_level = 1, SIZE(rs_rho_ext)
182 0 : CALL rs_grid_release(rs_rho_ext(igrid_level))
183 : END DO
184 0 : tot_rho_r_ext(1) = pw_integrate_function(rho_ext_r(1), isign=1)
185 0 : IF (my_rank == 0) THEN
186 0 : WRITE (*, FMT="(T3,A,T61,F20.10)") "ZMP| Total external charge: ", &
187 0 : tot_rho_r_ext(1)
188 : END IF
189 0 : DEALLOCATE (buffer, rs_rho_ext)
190 0 : CALL gid%sync()
191 : END ASSOCIATE
192 : END IF
193 :
194 10250 : CALL timestop(handle)
195 :
196 10250 : END SUBROUTINE external_read_density
197 :
198 : END MODULE qs_external_density
199 :
|
