Octopus
pcm.F90
Go to the documentation of this file.
1!! Copyright (C) 2014 Alain Delgado Gran, Carlo Andrea Rozzi, Stefano Corni, Gabriel Gil
2!!
3!! This program is free software; you can redistribute it and/or modify
4!! it under the terms of the GNU General Public License as published by
5!! the Free Software Foundation; either version 2, or (at your option)
6!! any later version.
7!!
8!! This program is distributed in the hope that it will be useful,
9!! but WITHOUT ANY WARRANTY; without even the implied warranty of
10!! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
11!! GNU General Public License for more details.
12!!
13!! You should have received a copy of the GNU General Public License
14!! along with this program; if not, write to the Free Software
15!! Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA
16!! 02110-1301, USA.
17!!
18
19#include "global.h"
20
21module pcm_oct_m
22 use box_minimum_oct_m
23 use debug_oct_m
24 use global_oct_m
25 use grid_oct_m
26 use, intrinsic :: iso_fortran_env
27 use index_oct_m
28 use io_oct_m
29 use ions_oct_m
30 use lalg_basic_oct_m
31 use math_oct_m
32 use messages_oct_m
33 use mesh_oct_m
34 use mesh_interpolation_oct_m
35 use mpi_oct_m
36 use namespace_oct_m
37 use par_vec_oct_m
38 use parser_oct_m
39 use pcm_eom_oct_m
40 use poisson_oct_m
41 use profiling_oct_m
42 use species_oct_m
43 use space_oct_m
44 use varinfo_oct_m
45
46 ! to output debug info
47 use unit_oct_m
48 use unit_system_oct_m
49 use io_function_oct_m
50 use mesh_function_oct_m
51
52 implicit none
53
54 private
55
56 public :: &
57 pcm_t, &
58 pcm_min_t, &
59 pcm_sphere_t, &
60 pcm_init, &
61 pcm_end, &
62 pcm_charges, &
63 pcm_get_vdw_radius, &
64 pcm_pot_rs, &
65 pcm_elect_energy, &
66 pcm_v_nuclei_cav, &
67 pcm_v_cav_li, &
68 pcm_update, &
69 pcm_calc_pot_rs, &
70 pcm_charge_density, &
71 pcm_dipole, &
72 pcm_eps, &
73 pcm_min_input_parsing_for_spectrum, &
74 pcm_td_eq, &
75 pcm_td_neq, &
76 pcm_td_eom
77
78
81 type pcm_sphere_t
82 private
83 real(real64) :: x !
84 real(real64) :: y
85 real(real64) :: z !
86 real(real64) :: r
87 end type pcm_sphere_t
88
90 integer, parameter :: PCM_DIM_SPACE = 3
91
92 type pcm_t
93 private
94 logical, public :: run_pcm
95 integer, public :: tdlevel
96 integer :: n_spheres
97 integer, public :: n_tesserae
98 type(pcm_sphere_t), allocatable :: spheres(:)
99 type(pcm_tessera_t), allocatable, public :: tess(:)
100 real(real64) :: scale_r
101 real(real64), allocatable :: matrix(:,:)
102 real(real64), allocatable :: matrix_d(:,:)
103 real(real64), allocatable :: matrix_lf(:,:)
104 real(real64), allocatable :: matrix_lf_d(:,:)
105 real(real64), allocatable, public :: q_e(:)
106 real(real64), allocatable :: q_n(:)
107 real(real64), allocatable, public :: q_e_in(:)
108 real(real64), allocatable :: rho_e(:)
109 real(real64), allocatable :: rho_n(:)
110 real(real64) :: qtot_e
111 real(real64) :: qtot_n
112 real(real64) :: qtot_e_in
113 real(real64) :: q_e_nominal
114 real(real64) :: q_n_nominal
115 logical :: renorm_charges
116 real(real64) :: q_tot_tol
117 real(real64) :: deltaQ_e
118 real(real64) :: deltaQ_n
119 real(real64), allocatable :: v_e(:)
120 real(real64), allocatable :: v_n(:)
121 real(real64), allocatable, public :: v_e_rs(:)
122 real(real64), allocatable, public :: v_n_rs(:)
123 real(real64), allocatable :: q_ext(:)
124 real(real64), allocatable :: q_ext_in(:)
125 real(real64), allocatable :: rho_ext(:)
126 real(real64) :: qtot_ext
127 real(real64) :: qtot_ext_in
128 real(real64), allocatable :: v_ext(:)
129 real(real64), allocatable, public :: v_ext_rs(:)
130 real(real64), allocatable :: q_kick(:)
131 real(real64), allocatable :: rho_kick(:)
132 real(real64) :: qtot_kick
133 real(real64), allocatable :: v_kick(:)
134 real(real64), allocatable, public :: v_kick_rs(:)
135 real(real64), public :: epsilon_0
136 real(real64), public :: epsilon_infty
137 integer, public :: which_eps
138 type(debye_param_t) :: deb
139 type(drude_param_t) :: drl
140 logical, public :: localf
141 logical, public :: solute
142 logical :: kick_is_present
143 ! !! (if localf is .false. this is irrelevant to PCM)
144 logical, public :: kick_like
145 integer :: initial_asc
146 real(real64) :: gaussian_width
147 integer :: info_unit
148 integer, public :: counter
149 character(len=80) :: input_cavity
150 integer :: update_iter
151 integer, public :: iter
152 integer :: calc_method
153 integer :: tess_nn
154 real(real64), public :: dt
155
156 ! We will allow a copy to the namespace and space here,
157 ! as this type will become an extension to the system_t class at some point
158 type(namespace_t), pointer :: namespace
159 type(space_t) :: space
160 end type pcm_t
161
162 type pcm_min_t
163 ! Components are public by default
164 logical :: run_pcm
165 logical :: localf
166 integer :: tdlevel
167 integer :: which_eps
168 type(debye_param_t) :: deb
169 type(drude_param_t) :: drl
170 end type pcm_min_t
171
172 real(real64), allocatable :: s_mat_act(:,:)
173 real(real64), allocatable :: d_mat_act(:,:)
174 real(real64), allocatable :: Sigma(:,:)
175 real(real64), allocatable :: Delta(:,:)
176
177 logical :: gamess_benchmark
178 real(real64), allocatable :: mat_gamess(:,:)
179
180 integer, parameter :: &
181 pcm_td_eq = 0, &
182 pcm_td_neq = 1, &
183 pcm_td_eom = 2
184
185 integer, parameter, public :: &
186 pcm_calc_direct = 1, &
187 pcm_calc_poisson = 2
188
189 integer, parameter :: &
190 pcm_vdw_optimized = 1, &
191 pcm_vdw_species = 2
192
193 integer, parameter :: n_tess_sphere = 60
194 ! !< cannot be changed without changing cav_gen subroutine
195
196contains
197
198 !-------------------------------------------------------------------------------------------------------
200 subroutine pcm_init(pcm, namespace, space, ions, grid, qtot, val_charge, external_potentials_present, kick_present)
201 type(pcm_t), intent(out) :: pcm
202 type(namespace_t), target, intent(in) :: namespace
203 class(space_t), intent(in) :: space
204 type(ions_t), intent(in) :: ions
205 type(grid_t), intent(in) :: grid
206 real(real64), intent(in) :: qtot
207 real(real64), intent(in) :: val_charge
208 logical, intent(in) :: external_potentials_present
209 logical, intent(in) :: kick_present
210
211 integer :: ia, ii, itess, jtess, pcm_vdw_type, subdivider
212 integer :: cav_unit_test, iunit, pcmmat_unit
213 integer :: pcmmat_gamess_unit, cav_gamess_unit
214 real(real64) :: min_distance
215
216 integer, parameter :: mxts = 10000
217
218 real(real64) :: default_value
219 real(real64) :: vdw_radius
220
221 type(pcm_tessera_t), allocatable :: dum2(:)
222
223 logical :: band
224 logical :: add_spheres_h
225 logical :: changed_default_nn
226
227 integer :: default_nn
228 real(real64) :: max_area
229
230 push_sub(pcm_init)
231
232 pcm%kick_is_present = kick_present
233
234 pcm%localf = .false.
235 pcm%iter = 0
236 pcm%update_iter = 1
237 pcm%kick_like = .false.
238
239 pcm%namespace => namespace
240 pcm%space = space
241
242 !%Variable PCMCalculation
243 !%Type logical
244 !%Default no
245 !%Section Hamiltonian::PCM
246 !%Description
247 !% If true, the calculation is performed accounting for solvation effects
248 !% by using the Integral Equation Formalism Polarizable Continuum Model IEF-PCM
249 !% formulated in real-space and real-time (<i>J. Chem. Phys.</i> <b>143</b>, 144111 (2015),
250 !% <i>Chem. Rev.</i> <b>105</b>, 2999 (2005), <i>J. Chem. Phys.</i> <b>139</b>, 024105 (2013)).
251 !% At the moment, this option is available only for <tt>TheoryLevel = DFT</tt>.
252 !% PCM is tested for CalculationMode = gs, while still experimental for other values (in particular, CalculationMode = td).
253 !%End
254 call parse_variable(namespace, 'PCMCalculation', .false., pcm%run_pcm)
255 if (pcm%run_pcm) then
256 call messages_print_with_emphasis(msg='PCM', namespace=namespace)
257 if (pcm%space%dim /= pcm_dim_space) then
258 message(1) = "PCM is only available for 3d calculations"
259 call messages_fatal(1, namespace=namespace)
260 end if
261 select type (box => grid%box)
262 type is (box_minimum_t)
263 class default
264 message(1) = "PCM is only available for BoxShape = minimum"
265 call messages_fatal(1, namespace=namespace)
266 end select
267 else
268 pop_sub(pcm_init)
269 return
270 end if
271
272 !%Variable PCMVdWRadii
273 !%Type integer
274 !%Default pcm_vdw_optimized
275 !%Section Hamiltonian::PCM
276 !%Description
277 !% This variable selects which van der Waals radius will be used to generate the solvent cavity.
278 !%Option pcm_vdw_optimized 1
279 !% Use the van der Waals radius optimized by Stefan Grimme in J. Comput. Chem. 27: 1787-1799, 2006,
280 !% except for C, N and O, reported in J. Chem. Phys. 120, 3893 (2004).
281 !%Option pcm_vdw_species 2
282 !% The vdW radii are set from the <tt>share/pseudopotentials/elements</tt> file. These values are obtained from
283 !% Alvarez S., Dalton Trans., 2013, 42, 8617-8636. Values can be changed in the <tt>Species</tt> block.
284 !%End
285 call parse_variable(namespace, 'PCMVdWRadii', pcm_vdw_optimized, pcm_vdw_type)
286 call messages_print_var_option("PCMVdWRadii", pcm_vdw_type, namespace=namespace)
287
288 select case (pcm_vdw_type)
289 case (pcm_vdw_optimized)
290 default_value = 1.2_real64
291 case (pcm_vdw_species)
292 default_value = m_one
293 end select
294
295 !%Variable PCMRadiusScaling
296 !%Type float
297 !%Section Hamiltonian::PCM
298 !%Description
299 !% Scales the radii of the spheres used to build the solute cavity surface.
300 !% The default value depends on the choice of <tt>PCMVdWRadii</tt>:
301 !% 1.2 for <tt>pcm_vdw_optimized</tt> and 1.0 for <tt>pcm_vdw_species</tt>.
302 !%End
303 call parse_variable(namespace, 'PCMRadiusScaling', default_value, pcm%scale_r)
304 call messages_print_var_value("PCMRadiusScaling", pcm%scale_r, namespace=namespace)
305
306 !%Variable PCMTDLevel
307 !%Type integer
308 !%Default eq
309 !%Section Hamiltonian::PCM
310 !%Description
311 !% When CalculationMode=td, PCMTDLevel it sets the way the time-depenendent solvent polarization is propagated.
312 !%Option eq 0
313 !% If PCMTDLevel=eq, the solvent is always in equilibrium with the solute or the external field, i.e.,
314 !% the solvent polarization follows instantaneously the changes in solute density or in the external field.
315 !% PCMTDLevel=neq and PCMTDLevel=eom are both nonequilibrium runs.
316 !%Option neq 1
317 !% If PCMTDLevel=neq, solvent polarization charges are splitted in two terms:
318 !% one that follows instantaneously the changes in the solute density or in the external field (dynamical polarization charges),
319 !% and another that lag behind in the evolution w.r.t. the solute density or the external field (inertial polarization charges).
320 !%Option eom 2
321 !% If PCMTDLevel=eom, solvent polarization charges evolves following an equation of motion, generalizing 'neq' propagation.
322 !% The equation of motion used here depends on the value of PCMEpsilonModel.
323 !%End
324 call parse_variable(namespace, 'PCMTDLevel', pcm_td_eq, pcm%tdlevel)
325 call messages_print_var_value("PCMTDLevel", pcm%tdlevel, namespace=namespace)
326
327 if (pcm%tdlevel /= pcm_td_eq .and. (.not. pcm%run_pcm)) then
328 call messages_write('Sorry, you have set PCMTDLevel /= eq, but PCMCalculation = no.')
329 call messages_new_line()
330 call messages_write('To spare you some time, Octopus will proceed as if PCMCalculation = yes.')
331 call messages_warning(namespace=namespace)
332 end if
333
334 !%Variable PCMStaticEpsilon
335 !%Type float
336 !%Default 1.0
337 !%Section Hamiltonian::PCM
338 !%Description
339 !% Static dielectric constant of the solvent (<math>\varepsilon_0</math>). 1.0 indicates gas phase.
340 !%End
341 call parse_variable(namespace, 'PCMStaticEpsilon', m_one, pcm%epsilon_0)
342 call messages_print_var_value("PCMStaticEpsilon", pcm%epsilon_0, namespace=namespace)
343
344 !%Variable PCMDynamicEpsilon
345 !%Type float
346 !%Default PCMStaticEpsilon
347 !%Section Hamiltonian::PCM
348 !%Description
349 !% High-frequency dielectric constant of the solvent (<math>\varepsilon_d</math>).
350 !% <math>\varepsilon_d=\varepsilon_0</math> indicate equilibrium with solvent.
351 !%End
352 call parse_variable(namespace, 'PCMDynamicEpsilon', pcm%epsilon_0, pcm%epsilon_infty)
353 call messages_print_var_value("PCMDynamicEpsilon", pcm%epsilon_infty, namespace=namespace)
354
355 !%Variable PCMEpsilonModel
356 !%Type integer
357 !%Default pcm_debye
358 !%Section Hamiltonian::PCM
359 !%Description
360 !% Define the dielectric function model.
361 !%Option pcm_debye 1
362 !% Debye model: <math>\varepsilon(\omega)=\varepsilon_d+\frac{\varepsilon_0-\varepsilon_d}{1-i\omega\tau}</math>
363 !%Option pcm_drude 2
364 !% Drude-Lorentz model: <math>\varepsilon(\omega)=1+\frac{A}{\omega_0^2-\omega^2+i\gamma\omega}</math>
365 !%End
366 call parse_variable(namespace, 'PCMEpsilonModel', pcm_debye_model, pcm%which_eps)
367 call messages_print_var_value("PCMEpsilonModel", pcm%which_eps, namespace=namespace)
368 if (.not. varinfo_valid_option('PCMEpsilonModel', pcm%which_eps)) call messages_input_error(namespace, 'PCMEpsilonModel')
369
370 if (pcm%tdlevel == pcm_td_eom .and. pcm%which_eps == pcm_drude_model) then
371 call messages_experimental("Drude-Lorentz EOM-PCM is experimental", namespace=namespace)
372 end if
373
374 if (pcm%which_eps /= pcm_debye_model .and. pcm%which_eps /= pcm_drude_model) then
375 call messages_write('Sorry, only Debye or Drude-Lorentz dielectric models are available.')
376 call messages_new_line()
377 call messages_write('To spare you some time, Octopus will proceed with the default choice (Debye).')
378 call messages_new_line()
379 call messages_write('You may change PCMEpsilonModel value for a Drude-Lorentz run.')
380 call messages_warning(namespace=namespace)
381 end if
382
383 if (pcm%tdlevel /= pcm_td_eq .and. pcm%which_eps == pcm_debye_model .and. is_close(pcm%epsilon_0, pcm%epsilon_infty)) then
384 call messages_write('Sorry, inertial/dynamic polarization splitting scheme for TD-PCM')
385 call messages_write('or Debye equation-of-motion TD-PCM')
386 call messages_new_line()
387 call messages_write('require both static and dynamic dielectric constants,')
388 call messages_write('and they must be different.')
389 call messages_new_line()
390 call messages_write('Octopus will run using TD-PCM version in equilibrium')
391 call messages_write('equilibrium with solvent at each time.')
392 call messages_warning(namespace=namespace)
393 pcm%tdlevel = pcm_td_eq
394 end if
395
396 !%Variable PCMEoMInitialCharges
397 !%Type integer
398 !%Default 0
399 !%Section Hamiltonian::PCM
400 !%Description
401 !% If =0 the propagation of the solvent polarization charges starts from internally generated initial charges
402 !% in equilibrium with the initial potential.
403 !% For Debye EOM-PCM, if >0 the propagation of the solvent polarization charges starts from initial charges from input file.
404 !% if =1, initial pol. charges due to solute electrons are read from input file.
405 !% else if =2, initial pol. charges due to external potential are read from input file.
406 !% else if =3, initial pol. charges due to solute electrons and external potential are read from input file.
407 !% Files should be located in pcm directory and are called ASC_e.dat and ASC_ext.dat, respectively.
408 !% The latter files are generated after any PCM run and contain the last values of the polarization charges.
409 !%End
410 call parse_variable(namespace, 'PCMEoMInitialCharges', 0, pcm%initial_asc)
411 call messages_print_var_value("PCMEoMInitialCharges", pcm%initial_asc, namespace=namespace)
412
413 if (pcm%initial_asc /= 0) then
414 call messages_experimental("The use of external initialization")
415 call messages_experimental("for the EOM-PCM charges is experimental", namespace=namespace)
416 if ((pcm%tdlevel /= pcm_td_eom .or. (pcm%tdlevel == pcm_td_eom .and. pcm%which_eps /= pcm_debye_model))) then
417 call messages_write('Sorry, initial polarization charges can only be read')
418 call messages_write('from input file for a Debye EOM-PCM run.')
419 call messages_new_line()
420 call messages_write('To spare you some time,')
421 call messages_write('Octopus will proceed as if PCMEoMInitialCharges = 0.')
422 call messages_warning(namespace=namespace)
423 pcm%initial_asc = 0
424 end if
425 end if
426
428 pcm%deb%eps_0 = pcm%epsilon_0
429 pcm%deb%eps_d = pcm%epsilon_infty
430
432 call parse_variable(namespace, 'TDTimeStep', m_zero, pcm%dt, unit = units_inp%time)
433
434 !%Variable PCMDebyeRelaxTime
435 !%Type float
436 !%Default 0.0
437 !%Section Hamiltonian::PCM
438 !%Description
439 !% Relaxation time of the solvent within Debye model (<math>\tau</math>). Recall Debye dieletric function:
440 !% <math>\varepsilon(\omega)=\varepsilon_d+\frac{\varepsilon_0-\varepsilon_d}{1-i\omega\tau}</math>
441 !%End
442 call parse_variable(namespace, 'PCMDebyeRelaxTime', m_zero, pcm%deb%tau)
443 call messages_print_var_value("PCMDebyeRelaxTime", pcm%deb%tau, namespace=namespace)
444
445 if (pcm%tdlevel == pcm_td_eom .and. pcm%which_eps == pcm_debye_model .and. &
446 (abs(pcm%deb%tau) <= m_epsilon .or. is_close(pcm%deb%eps_0, pcm%deb%eps_d))) then
447 call messages_write('Sorry, you have set PCMTDLevel = eom,')
448 call messages_write('but you have not included all required Debye model parameters.')
449 call messages_new_line()
450 call messages_write('You need PCMEpsilonStatic, PCMEpsilonDynamic')
451 call messages_write('and PCMDebyeRelaxTime for an EOM TD-PCM run.')
452 call messages_new_line()
453 call messages_write('Octopus will run using TD-PCM version in equilibrium')
454 call messages_write('equilibrium with solvent at each time.')
455 call messages_warning(namespace=namespace)
456 pcm%tdlevel = pcm_td_eq
457 end if
458
459 if (is_close(pcm%epsilon_0, m_one)) then
460 if (pcm%tdlevel == pcm_td_eom .and. pcm%which_eps == pcm_drude_model) then
461 message(1) = "PCMEpsilonStatic = 1 is incompatible with a Drude-Lorentz EOM-PCM run."
462 call messages_fatal(1, namespace=namespace)
463 end if
464 else
465 !%Variable PCMDrudeLOmega
466 !%Type float
467 !%Default <math>\sqrt{1/(\varepsilon_0-1)}</math>
468 !%Section Hamiltonian::PCM
469 !%Description
470 !% Resonance frequency of the solvent within Drude-Lorentz model (<math>\omega_0</math>).
471 !% Recall Drude-Lorentz dielectric function: <math>\varepsilon(\omega)=1+\frac{A}{\omega_0^2-\omega^2+i\gamma\omega}</math>
472 !% Default values of <math>\omega_0</math> guarantee to recover static dielectric constant.
473 !%End
474 call parse_variable(namespace, 'PCMDrudeLOmega', sqrt(m_one/(pcm%epsilon_0 - m_one)), pcm%drl%w0)
475 call messages_print_var_value("PCMDrudeLOmega", pcm%drl%w0, namespace=namespace)
476 end if
477
478 if (pcm%tdlevel == pcm_td_eom .and. pcm%which_eps == pcm_drude_model .and. abs(pcm%drl%w0) <= m_epsilon) then
479 call messages_write('Sorry, you have set PCMDrudeLOmega =')
480 call messages_write('but this is incompatible with a Drude-Lorentz EOM-PCM run.')
481 call messages_new_line()
482 if (.not. is_close(pcm%epsilon_0, m_one)) then
483 call messages_write('Octopus will run using the default value of PCMDrudeLOmega.')
484 call messages_warning(namespace=namespace)
485 pcm%drl%w0 = sqrt(m_one/(pcm%epsilon_0 - m_one))
486 else
487 message(1) = "PCMEpsilonStatic = 1 is incompatible with a Drude-Lorentz EOM-PCM run."
488 call messages_fatal(1, namespace=namespace)
489 end if
490 end if
491
492 !%Variable PCMDrudeLDamping
493 !%Type float
494 !%Default 0.0
495 !%Section Hamiltonian::PCM
496 !%Description
497 !% Damping factor of the solvent charges oscillations within Drude-Lorentz model (<math>\gamma</math>).
498 !% Recall Drude-Lorentz dielectric function: <math>\varepsilon(\omega)=1+\frac{A}{\omega_0^2-\omega^2+i\gamma\omega}</math>
499 !%End
500 call parse_variable(namespace, 'PCMDrudeLDamping', m_zero, pcm%drl%gm)
501 call messages_print_var_value("PCMDrudeLDamping", pcm%drl%gm, namespace=namespace)
502
503
504 pcm%drl%aa = (pcm%epsilon_0 - m_one)*pcm%drl%w0**2
505
506 !%Variable PCMLocalField
507 !%Type logical
508 !%Default no
509 !%Section Hamiltonian::PCM
510 !%Description
511 !% This variable is a flag for including local field effects when an external field is applied. The total field interacting with
512 !% the molecule (also known as cavity field) is not the bare field in the solvent (the so-called Maxwell field), but it also
513 !% include a contribution due to the polarization of the solvent. The latter is calculated here within the PCM framework.
514 !% See [G. Gil, et al., J. Chem. Theory Comput., 2019, 15 (4), pp 2306–2319].
515 !%End
516 call parse_variable(namespace, 'PCMLocalField', .false., pcm%localf)
517 call messages_print_var_value("PCMLocalField", pcm%localf, namespace=namespace)
518
519 if (pcm%localf .and. ((.not. external_potentials_present) .and. (.not. pcm%kick_is_present))) then
520 message(1) = "Sorry, you have set PCMLocalField = yes, but you have not included any external potentials."
521 call messages_fatal(1, namespace=namespace)
522 end if
523
524 !%Variable PCMSolute
525 !%Type logical
526 !%Default yes
527 !%Section Hamiltonian::PCM
528 !%Description
529 !% This variable is a flag for including polarization effects of the solvent due to the solute.
530 !% (Useful for analysis) When external fields are applied, turning off the solvent-molecule interaction (PCMSolute=no) and
531 !% activating the solvent polarization due to the applied field (PCMLocalField=yes) allows to include only local field effects.
532 !%End
533 call parse_variable(namespace, 'PCMSolute', .true., pcm%solute)
534 call messages_print_var_value("PCMSolute", pcm%solute, namespace=namespace)
535
536 if (pcm%run_pcm .and. (.not. pcm%solute)) then
537 call messages_write('N.B. This PCM run do not consider the polarization effects due to the solute.')
538 call messages_warning(namespace=namespace)
539 if (.not. pcm%localf) then
540 message(1) = "You have activated a PCM run without polarization effects. Octopus will halt."
541 call messages_fatal(1, namespace=namespace)
542 end if
543 end if
544
545 !%Variable PCMKick
546 !%Type logical
547 !%Default no
548 !%Section Hamiltonian::PCM
549 !%Description
550 !% This variable controls the effect the kick has on the polarization of the solvent.
551 !% If .true. ONLY the FAST degrees-of-freedom of the solvent follow the kick. The potential due to polarization charges behaves
552 !% as another kick, i.e., it is a delta-perturbation.
553 !% If .false. ALL degrees-of-freedom of the solvent follow the kick. The potential due to polarization charges evolves
554 !% following an equation of motion. When Debye dielectric model is used, just a part of the potential behaves as another kick.
555 !%End
556 call parse_variable(namespace, 'PCMKick', .false., pcm%kick_like)
557 call messages_print_var_value("PCMKick", pcm%kick_like, namespace=namespace)
558
559 if (pcm%kick_like .and. (.not. pcm%run_pcm)) then
560 message(1) = "PCMKick option can only be activated when PCMCalculation = yes. Octopus will halt."
561 call messages_fatal(1, namespace=namespace)
562 end if
563
564 if (pcm%kick_like .and. (.not. pcm%localf)) then
565 message(1) = "PCMKick option can only be activated when a PCMLocalField = yes. Octopus will halt."
566 call messages_fatal(1, namespace=namespace)
567 end if
568
569 if (pcm%kick_like .and. (.not. pcm%kick_is_present)) then
570 message(1) = "Sorry, you have set PCMKick = yes, but you have not included any kick."
571 call messages_fatal(1, namespace=namespace)
572 end if
573
574 if (pcm%kick_is_present .and. pcm%run_pcm .and. (.not. pcm%localf)) then
575 message(1) = "You have set up a PCM calculation with a kick without local field effects."
576 message(2) = "Please, reconsider if you want the kick to be relevant for the PCM run."
577 call messages_warning(2, namespace=namespace)
578 end if
579
580 !%Variable PCMUpdateIter
581 !%Type integer
582 !%Default 1
583 !%Section Hamiltonian::PCM
584 !%Description
585 !% Defines how often the PCM potential is updated during time propagation.
586 !%End
587 call parse_variable(namespace, 'PCMUpdateIter', 1, pcm%update_iter)
588 call messages_print_var_value("PCMUpdateIter", pcm%update_iter, namespace=namespace)
589
590 !%Variable PCMGamessBenchmark
591 !%Type logical
592 !%Default .false.
593 !%Section Hamiltonian::PCM
594 !%Description
595 !% If PCMGamessBenchmark is set to "yes", the pcm_matrix is also written in a Gamess format.
596 !% for benchamarking purposes.
597 !%End
598 call parse_variable(namespace, 'PCMGamessBenchmark', .false., gamess_benchmark)
599
600 !%Variable PCMRenormCharges
601 !%Type logical
602 !%Default .false.
603 !%Section Hamiltonian::PCM
604 !%Description
605 !% If .true. renormalization of the polarization charges is performed to enforce fulfillment
606 !% of the Gauss law, <math>\sum_i q_i^{e/n} = -[(\epsilon-1)/\epsilon] Q_M^{e/n}</math> where
607 !% <math>q_i^{e/n}</math> are the polarization charges induced by the electrons/nuclei of the molecule
608 !% and <math>Q_M^{e/n}</math> is the nominal electronic/nuclear charge of the system. This can be needed
609 !% to treat molecules in weakly polar solvents.
610 !%End
611 call parse_variable(namespace, 'PCMRenormCharges', .false., pcm%renorm_charges)
612
613 !%Variable PCMQtotTol
614 !%Type float
615 !%Default 0.5
616 !%Section Hamiltonian::PCM
617 !%Description
618 !% If <tt>PCMRenormCharges=.true.</tt> and <math>\delta Q = |[\sum_i q_i| - ((\epsilon-1)/\epsilon)*|Q_M]|>PCMQtotTol</math>
619 !% the polarization charges will be normalized as
620 !% <math>q_i^\prime=q_i + signfunction(e, n, \delta Q) (q_i/q_{tot})*\delta Q</math>
621 !% with <math>q_{tot} = \sum_i q_i</math>. For values of <math>\delta Q > 0.5</math>
622 !% (printed by the code in the file pcm/pcm_info.out) even, if polarization charges are renormalized,
623 !% the calculated results might be inaccurate or erroneous.
624 !%End
625 call parse_variable(namespace, 'PCMQtotTol', m_half, pcm%q_tot_tol)
626
627 if (pcm%renorm_charges) then
628 message(1) = "Info: Polarization charges will be renormalized"
629 message(2) = " if |Q_tot_PCM - Q_M| > PCMQtotTol"
630 call messages_info(2, namespace=namespace)
631 end if
632
633 !%Variable PCMSmearingFactor
634 !%Type float
635 !%Default 1.0
636 !%Section Hamiltonian::PCM
637 !%Description
638 !% Parameter used to control the width (area of each tessera) of the Gaussians used to distribute
639 !% the polarization charges on each tessera (arXiv:1507.05471). If set to zero, the solvent
640 !% reaction potential in real-space is defined by using point charges.
641 !%End
642 call parse_variable(namespace, 'PCMSmearingFactor', m_one, pcm%gaussian_width)
643 call messages_print_var_value("PCMSmearingFactor", pcm%gaussian_width, namespace=namespace)
644
645 if (abs(pcm%gaussian_width) <= m_epsilon) then
646 message(1) = "Info: PCM potential will be defined in terms of polarization point charges"
647 call messages_info(1, namespace=namespace)
648 else
649 message(1) = "Info: PCM potential is regularized to avoid Coulomb singularity"
650 call messages_info(1, namespace=namespace)
651 end if
652
653 call io_mkdir('pcm', namespace)
654
655 !%Variable PCMCavity
656 !%Type string
657 !%Section Hamiltonian::PCM
658 !%Description
659 !% Name of the file containing the geometry of the cavity hosting the solute molecule.
660 !% The data must be in atomic units and the file must contain the following information sequentially:
661 !% T < Number of tesserae
662 !% s_x(1:T) < coordinates x of the tesserae
663 !% s_y(1:T) < coordinates y of the tesserae
664 !% s_z(1:T) < coordinates z of the tesserae
665 !% A(1:T) < areas of the tesserae
666 !% R_sph(1:T) < Radii of the spheres to which the tesserae belong
667 !% normal(1:T,1:3) < Outgoing unitary vectors at the tesserae surfaces
668 !%End
669 call parse_variable(namespace, 'PCMCavity', '', pcm%input_cavity)
670
671 if (pcm%input_cavity == '') then
672
673 !%Variable PCMSpheresOnH
674 !%Type logical
675 !%Default no
676 !%Section Hamiltonian::PCM
677 !%Description
678 !% If true, spheres centered at the Hydrogens atoms are included to build the solute cavity surface.
679 !%End
680 call parse_variable(namespace, 'PCMSpheresOnH', .false., add_spheres_h)
681
682 pcm%n_spheres = 0
683 band = .false.
684 do ia = 1, ions%natoms
685 if ((.not. add_spheres_h) .and. ions%atom(ia)%label == 'H') cycle
686 pcm%n_spheres = pcm%n_spheres + 1 !counting the number of species different from Hydrogen
687 end do
688
689 safe_allocate(pcm%spheres(1:pcm%n_spheres))
690 pcm%spheres(:)%x = m_zero
691 pcm%spheres(:)%y = m_zero
692 pcm%spheres(:)%z = m_zero
693 pcm%spheres(:)%r = m_zero
694
695 pcm%n_spheres = 0
696 do ia = 1, ions%natoms
697 if ((.not. add_spheres_h) .and. ions%atom(ia)%label == 'H') cycle
698 pcm%n_spheres = pcm%n_spheres + 1
699
701 pcm%spheres(pcm%n_spheres)%x = ions%pos(1, ia)
702 pcm%spheres(pcm%n_spheres)%y = ions%pos(2, ia)
703 pcm%spheres(pcm%n_spheres)%z = ions%pos(3, ia)
704
705 vdw_radius = pcm_get_vdw_radius(ions%atom(ia)%species, pcm_vdw_type, namespace)
706 pcm%spheres(pcm%n_spheres)%r = vdw_radius*pcm%scale_r
707 end do
708
709 if (mpi_world%is_root()) then
710 pcm%info_unit = io_open(pcm_dir//'pcm_info.out', pcm%namespace, action='write')
711
712 write(pcm%info_unit, '(A35)') '# Configuration: Molecule + Solvent'
713 write(pcm%info_unit, '(A35)') '# ---------------------------------'
714 write(pcm%info_unit, '(A21,F12.3)') '# Epsilon(Solvent) = ', pcm%epsilon_0
715 write(pcm%info_unit, '(A1)')'#'
716 write(pcm%info_unit, '(A35,I4)') '# Number of interlocking spheres = ', pcm%n_spheres
717 write(pcm%info_unit, '(A1)')'#'
718
719 write(pcm%info_unit, '(A8,3X,A7,8X,A26,20X,A10)') '# SPHERE', 'ELEMENT', 'CENTER (X,Y,Z) (A)', 'RADIUS (A)'
720 write(pcm%info_unit, '(A8,3X,A7,4X,A43,7X,A10)') '# ------', '-------', &
721 '-------------------------------------------', '----------'
722 end if
723
724 pcm%n_spheres = 0
725 do ia = 1, ions%natoms
726 if ((.not. add_spheres_h) .and. ions%atom(ia)%label == 'H') cycle
727 pcm%n_spheres = pcm%n_spheres + 1
728 if (mpi_world%is_root()) then
729 write(pcm%info_unit,'(A1,2X,I3,7X,A2,3X,F14.8,2X,F14.8,2X,F14.8,4X,F14.8)')'#', pcm%n_spheres, &
730 ions%atom(ia)%label, &
731 ions%pos(1, ia)*p_a_b, &
732 ions%pos(2, ia)*p_a_b, &
733 ions%pos(3, ia)*p_a_b, &
734 pcm%spheres(pcm%n_spheres)%r*p_a_b
735 end if
736 end do
737
738 !%Variable PCMTessSubdivider
739 !%Type integer
740 !%Default 1
741 !%Section Hamiltonian::PCM
742 !%Description
743 !% Allows to subdivide further each tessera refining the discretization of the cavity tesselation.
744 !% Can take only two values, 1 or 4. 1 corresponds to 60 tesserae per sphere, while 4 corresponds to 240 tesserae per sphere.
745 !%End
746 call parse_variable(namespace, 'PCMTessSubdivider', 1, subdivider)
747
748 safe_allocate(dum2(1:subdivider*n_tess_sphere*pcm%n_spheres))
749
750 !%Variable PCMTessMinDistance
751 !%Type float
752 !%Default 0.1
753 !%Section Hamiltonian::PCM
754 !%Description
755 !% Minimum distance between tesserae.
756 !% Any two tesserae having smaller distance in the starting tesselation will be merged together.
757 !%End
758 call parse_variable(namespace, 'PCMTessMinDistance', 0.1_real64, min_distance)
759 call messages_print_var_value("PCMTessMinDistance", min_distance, namespace=namespace)
760
762 call cav_gen(subdivider, min_distance, pcm%n_spheres, pcm%spheres, pcm%n_tesserae, dum2, pcm%info_unit)
763
764 safe_allocate(pcm%tess(1:pcm%n_tesserae))
765
766 do ia=1, pcm%n_tesserae
767 pcm%tess(ia)%point = m_zero
768 pcm%tess(ia)%area = m_zero
769 pcm%tess(ia)%r_sphere = m_zero
770 pcm%tess(ia)%normal = m_zero
771 end do
772
773 pcm%tess = dum2(1:pcm%n_tesserae)
774
775 safe_deallocate_a(dum2)
776
777 message(1) = "Info: van der Waals surface has been calculated"
778 call messages_info(1, namespace=namespace)
779
780 else
781
783 iunit = io_open(trim(pcm%input_cavity), pcm%namespace, action='read', status='old')
784 read(iunit,*) pcm%n_tesserae
785
786 if (pcm%n_tesserae > mxts) then
787 write(message(1),'(a,I5,a,I5)') "total number of tesserae", pcm%n_tesserae, ">", mxts
788 call messages_warning(1, namespace=namespace)
789 end if
790
791 safe_allocate(pcm%tess(1:pcm%n_tesserae))
792
793 do ia = 1, pcm%n_tesserae
794 pcm%tess(ia)%point = m_zero
795 pcm%tess(ia)%area = m_zero
796 pcm%tess(ia)%r_sphere = m_zero
797 pcm%tess(ia)%normal = m_zero
798 end do
799
800 do ia = 1, pcm%n_tesserae
801 read(iunit,*) pcm%tess(ia)%point(1)
802 end do
803
804 do ia = 1, pcm%n_tesserae
805 read(iunit,*) pcm%tess(ia)%point(2)
806 end do
807
808 do ia = 1, pcm%n_tesserae
809 read(iunit,*) pcm%tess(ia)%point(3)
810 end do
811
812 do ia = 1, pcm%n_tesserae
813 read(iunit,*) pcm%tess(ia)%area
814 end do
815
816 do ia = 1, pcm%n_tesserae
817 read(iunit,*) pcm%tess(ia)%r_sphere
818 end do
819
820 do ia = 1, pcm%n_tesserae
821 read(iunit,*) pcm%tess(ia)%normal
822 end do
823
824 call io_close(iunit)
825 message(1) = "Info: van der Waals surface has been read from " // trim(pcm%input_cavity)
826 call messages_info(1, namespace=namespace)
827 end if
828
829 if (mpi_world%is_root()) then
830 cav_unit_test = io_open(pcm_dir//'cavity_mol.xyz', pcm%namespace, action='write')
831
832 write (cav_unit_test,'(2X,I4)') pcm%n_tesserae + ions%natoms
833 write (cav_unit_test,'(2X)')
834
835 do ia = 1, pcm%n_tesserae
836 write(cav_unit_test,'(2X,A2,3X,4f15.8,3X,4f15.8,3X,4f15.8)') 'H', pcm%tess(ia)%point*p_a_b
837 end do
838
839 do ia = 1, ions%natoms
840 write(cav_unit_test,'(2X,A2,3X,4f15.8,3X,4f15.8,3X,4f15.8)') ions%atom(ia)%label, &
841 ions%pos(:, ia)*p_a_b
842 end do
843
844 call io_close(cav_unit_test)
845
846 write(pcm%info_unit,'(A1)')'#'
847 if (pcm%localf) then
848 write(pcm%info_unit,'(A1,4X,A4,14X,A4,21X,A4,21X,A4,21X,A4,21X,A7,18X,A7,18X,A8,17X,A5,20X,A8,17X,A5,20X,A8,17X,A5)') &
849 '#','iter', 'E_ee', 'E_en', 'E_nn', 'E_ne', 'E_e_ext', 'E_n_ext', 'E_M-solv', &
850 'Q_pol^e', 'deltaQ^e', 'Q_pol^n', 'deltaQ^n', 'Q_pol^ext'
851 else
852 write(pcm%info_unit,'(A1,4X,A4,14X,A4,21X,A4,21X,A4,21X,A4,21X,A8,17X,A5,20X,A8,17X,A5,20X, A8)') &
853 '#','iter', 'E_ee', 'E_en', 'E_nn', 'E_ne', 'E_M-solv', 'Q_pol^e', 'deltaQ^e', 'Q_pol^n', 'deltaQ^n'
854 end if
855 end if
856 pcm%counter = 0
857
859 if (gamess_benchmark .and. mpi_world%is_root()) then
860 cav_gamess_unit = io_open(pcm_dir//'geom_cavity_gamess.out', pcm%namespace, action='write')
861
862 write(cav_gamess_unit,*) pcm%n_tesserae
863
864 do ia = 1, pcm%n_tesserae
865 write(cav_gamess_unit,*) pcm%tess(ia)%point(1)
866 end do
867
868 do ia = 1, pcm%n_tesserae
869 write(cav_gamess_unit,*) pcm%tess(ia)%point(2)
870 end do
871
872 do ia = 1, pcm%n_tesserae
873 write(cav_gamess_unit,*) pcm%tess(ia)%point(3)
874 end do
875
876 do ia = 1, pcm%n_tesserae
877 write(cav_gamess_unit,*) pcm%tess(ia)%area
878 end do
879
880 do ia = 1, pcm%n_tesserae
881 write(cav_gamess_unit,*) pcm%tess(ia)%r_sphere
882 end do
883
884 do ia = 1, pcm%n_tesserae
885 write(cav_gamess_unit,*) pcm%tess(ia)%normal
886 end do
887
888 call io_close(cav_gamess_unit)
889 end if
890
892 if (gamess_benchmark) then
893 safe_allocate(mat_gamess(1:pcm%n_tesserae, 1:pcm%n_tesserae))
894 mat_gamess = m_zero
895 end if
896
897 if (.not. is_close(pcm%epsilon_infty, pcm%epsilon_0)) then
898
899 safe_allocate(pcm%matrix_d(1:pcm%n_tesserae, 1:pcm%n_tesserae))
900 pcm%matrix_d = m_zero
901
902 call pcm_matrix(pcm%epsilon_infty, pcm%tess, pcm%n_tesserae, pcm%matrix_d)
903
904 if (gamess_benchmark .and. mpi_world%is_root()) then
905 pcmmat_gamess_unit = io_open(pcm_dir//'pcm_matrix_gamess_dyn.out', pcm%namespace, action='write')
906
907 do jtess = 1, pcm%n_tesserae
908 do itess = 1, pcm%n_tesserae
909 write(pcmmat_gamess_unit,*) mat_gamess(itess,jtess)
910 end do
911 end do
912
913 call io_close(pcmmat_gamess_unit)
914 mat_gamess = m_zero
915 end if
916
917 if (pcm%localf) then
918
919 safe_allocate(pcm%matrix_lf_d(1:pcm%n_tesserae, 1:pcm%n_tesserae))
920 pcm%matrix_lf_d = m_zero
921
922 call pcm_matrix(pcm%epsilon_infty, pcm%tess, pcm%n_tesserae, pcm%matrix_lf_d, .true.)
923
924 if (mpi_world%is_root()) then
925 ! only to benchmark
926 pcmmat_unit = io_open(pcm_dir//'pcm_matrix_dynamic_lf.out', pcm%namespace, action='write')
927 do jtess = 1, pcm%n_tesserae
928 do itess = 1, pcm%n_tesserae
929 write(pcmmat_unit,*) pcm%matrix_lf_d(itess,jtess)
930 end do
931 end do
932 call io_close(pcmmat_unit)
933 end if
934
935 end if
936
937 end if
938
939 safe_allocate(pcm%matrix(1:pcm%n_tesserae, 1:pcm%n_tesserae))
940 pcm%matrix = m_zero
941
942 call pcm_matrix(pcm%epsilon_0, pcm%tess, pcm%n_tesserae, pcm%matrix)
943
944 if (mpi_world%is_root()) then
945 pcmmat_unit = io_open(pcm_dir//'pcm_matrix.out', pcm%namespace, action='write')
946 if (gamess_benchmark) pcmmat_gamess_unit = io_open(pcm_dir//'pcm_matrix_gamess.out', &
947 pcm%namespace, action='write')
948
949 do jtess = 1, pcm%n_tesserae
950 do itess = 1, pcm%n_tesserae
951 write(pcmmat_unit,*) pcm%matrix(itess,jtess)
952 if (gamess_benchmark) write(pcmmat_gamess_unit,*) mat_gamess(itess,jtess)
953 end do
954 end do
955 call io_close(pcmmat_unit)
956 if (gamess_benchmark) call io_close(pcmmat_gamess_unit)
957 end if
958 call mpi_world%barrier()
959
960 if (gamess_benchmark) then
961 safe_deallocate_a(mat_gamess)
962 end if
963
964 if (pcm%localf) then
965
966 safe_allocate(pcm%matrix_lf(1:pcm%n_tesserae, 1:pcm%n_tesserae))
967 pcm%matrix_lf = m_zero
968
969 call pcm_matrix(pcm%epsilon_0, pcm%tess, pcm%n_tesserae, pcm%matrix_lf, .true.)
970
971 if (mpi_world%is_root()) then
972 ! only to benchmark
973 pcmmat_unit = io_open(pcm_dir//'pcm_matrix_static_lf.out', pcm%namespace, action='write')
974 do jtess = 1, pcm%n_tesserae
975 do itess = 1, pcm%n_tesserae
976 write(pcmmat_unit,*) pcm%matrix_lf(itess,jtess)
977 end do
978 end do
979 call io_close(pcmmat_unit)
980 end if
981
982 end if
983
984 message(1) = "Info: PCM response matrices has been evaluated"
985 call messages_info(1, namespace=namespace)
986
987 !%Variable PCMCalcMethod
988 !%Type integer
989 !%Default pcm_direct
990 !%Section Hamiltonian::PCM
991 !%Description
992 !% Defines the method to be used to obtain the PCM potential.
993 !%Option pcm_direct 1
994 !% Direct sum of the potential generated by the polarization charges regularized
995 !% with a Gaussian smearing [A. Delgado, et al., J Chem Phys 143, 144111 (2015)].
996 !%Option pcm_poisson 2
997 !% Solving the Poisson equation for the polarization charge density.
998 !%End
999 call parse_variable(namespace, 'PCMCalcMethod', pcm_calc_direct, pcm%calc_method)
1000 call messages_print_var_option("PCMCalcMethod", pcm%calc_method, namespace=namespace)
1001
1002
1003 if (pcm%calc_method == pcm_calc_poisson) then
1004
1005 max_area = m_epsilon
1006 do ia = 1, pcm%n_tesserae
1007 if (pcm%tess(ia)%area > max_area) max_area = pcm%tess(ia)%area
1008 end do
1009
1010 !default is as many neighbor to contain 1 gaussian width
1011 default_nn = int(max_area*pcm%gaussian_width/minval(grid%spacing(1:pcm%space%dim)))
1012
1013 changed_default_nn = .false.
1014
1015 do ii = default_nn, 1, -1
1016 pcm%tess_nn = ii
1017 if (pcm_nn_in_mesh(pcm,grid)) then
1018 exit
1019 else
1020 changed_default_nn = .true.
1021 end if
1022 end do
1023 if (changed_default_nn) then
1024 call messages_write('PCM nearest neighbors have been reduced from ')
1025 call messages_write(default_nn)
1026 call messages_write(' to ')
1027 call messages_write(pcm%tess_nn)
1028 call messages_new_line()
1029 call messages_write('in order to fit them into the mesh.')
1030 call messages_new_line()
1031 call messages_write('This may produce unexpected results. ')
1032 call messages_warning(namespace=namespace)
1033 end if
1034 default_nn = pcm%tess_nn
1035
1036 !%Variable PCMChargeSmearNN
1037 !%Type integer
1038 !%Default 2 * max_area * PCMSmearingFactor
1039 !%Section Hamiltonian::PCM
1040 !%Description
1041 !% Defines the number of nearest neighbor mesh-points to be taken around each
1042 !% cavity tessera in order to smear the charge when PCMCalcMethod = pcm_poisson.
1043 !% Setting PCMChargeSmearNN = 1 means first nearest neighbors, PCMChargeSmearNN = 2
1044 !% second nearest neighbors, and so on.
1045 !% The default value is such that the neighbor mesh contains points in a radius
1046 !% equal to the width used for the gaussian smearing.
1047 !%End
1048 call parse_variable(namespace, 'PCMChargeSmearNN', default_nn, pcm%tess_nn)
1049 call messages_print_var_value("PCMChargeSmearNN", pcm%tess_nn, namespace=namespace)
1050
1051 call pcm_poisson_sanity_check(pcm, grid)
1052
1053 end if
1054
1055 if (pcm%run_pcm) call messages_print_with_emphasis(namespace=namespace)
1056
1057 if (pcm%calc_method == pcm_calc_poisson) then
1058 safe_allocate(pcm%rho_n(1:grid%np_part))
1059 safe_allocate(pcm%rho_e(1:grid%np_part))
1060 if (pcm%localf) then
1061 safe_allocate(pcm%rho_ext(1:grid%np_part))
1062 if (pcm%kick_is_present) then
1063 safe_allocate(pcm%rho_kick(1:grid%np_part))
1064 end if
1065 end if
1066 end if
1067
1068
1069 safe_allocate(pcm%v_n(1:pcm%n_tesserae))
1070 safe_allocate(pcm%q_n(1:pcm%n_tesserae))
1071 safe_allocate(pcm%v_n_rs(1:grid%np))
1072 pcm%v_n = m_zero
1073 pcm%q_n = m_zero
1074 pcm%v_n_rs = m_zero
1075
1076 safe_allocate(pcm%v_e(1:pcm%n_tesserae))
1077 safe_allocate(pcm%q_e(1:pcm%n_tesserae))
1078 safe_allocate(pcm%v_e_rs(1:grid%np))
1079 pcm%v_e = m_zero
1080 pcm%q_e = m_zero
1081 pcm%v_e_rs = m_zero
1082
1083 safe_allocate(pcm%q_e_in(1:pcm%n_tesserae))
1084 pcm%q_e_in = m_zero
1085
1086
1087 if (pcm%localf) then
1088 safe_allocate(pcm%v_ext(1:pcm%n_tesserae))
1089 safe_allocate(pcm%q_ext(1:pcm%n_tesserae))
1090 safe_allocate(pcm%v_ext_rs(1:grid%np))
1091 pcm%v_ext = m_zero
1092 pcm%q_ext = m_zero
1093 pcm%v_ext_rs = m_zero
1094
1095 safe_allocate(pcm%q_ext_in(1:pcm%n_tesserae))
1096 pcm%q_ext_in = m_zero
1097
1098 if (pcm%kick_is_present) then
1099 safe_allocate(pcm%v_kick(1:pcm%n_tesserae))
1100 safe_allocate(pcm%q_kick(1:pcm%n_tesserae))
1101 safe_allocate(pcm%v_kick_rs(1:grid%np))
1102 pcm%v_ext = m_zero
1103 pcm%q_ext = m_zero
1104 pcm%v_ext_rs = m_zero
1105 end if
1106
1107 end if
1108
1109 pcm%q_e_nominal = qtot
1110 pcm%q_n_nominal = val_charge
1111 pcm%deltaQ_e = m_zero
1112 pcm%deltaQ_n = m_zero
1113 pcm%qtot_e = m_zero
1114 pcm%qtot_n = m_zero
1115 pcm%qtot_ext = m_zero
1116
1117 pop_sub(pcm_init)
1118 end subroutine pcm_init
1119
1120 ! -----------------------------------------------------------------------------
1121 subroutine pcm_calc_pot_rs(pcm, mesh, psolver, ions, v_h, v_ext, kick, time_present, kick_time)
1122 type(pcm_t), intent(inout) :: pcm
1123 class(mesh_t), intent(in) :: mesh
1124 type(poisson_t), intent(in) :: psolver
1125 type(ions_t), optional, intent(in) :: ions
1126 real(real64), optional, intent(in) :: v_h(:)
1127 real(real64), optional, intent(in) :: v_ext(:)
1128 real(real64), optional, intent(in) :: kick(:)
1129 logical, optional, intent(in) :: time_present
1130 logical, optional, intent(in) :: kick_time
1131
1132 integer, save :: calc
1133
1134 logical, save :: input_asc_e
1135 logical, save :: input_asc_ext
1136
1137 ! for debuggin - it should be cleaned up
1138 !character*5 :: iteration
1139
1140 logical, save :: not_yet_called = .false.
1141 logical, save :: is_time_for_kick = .false.
1142 logical, save :: after_kick = .false.
1143
1144 logical, save :: td_calc_mode = .false.
1145
1146 integer, save :: asc_vs_t_unit_e
1147
1148 character(len=23), save :: asc_vs_t_unit_format
1149 character(len=16), save :: asc_vs_t_unit_format_tail
1150
1151 integer, save :: ia
1152
1153 push_sub(pcm_calc_pot_rs)
1154
1155 assert(present(v_h) .or. present(ions) .or. present(v_ext) .or. present(kick))
1156
1157 if (present(v_h)) calc = pcm_electrons
1158 if (present(ions)) calc = pcm_nuclei
1159 if (present(v_ext) .and. (.not. present(kick))) calc = pcm_external_potential
1160 if ((.not. present(v_ext)) .and. present(kick)) calc = pcm_kick
1161 if (present(v_ext) .and. present(kick)) calc = pcm_external_plus_kick
1162
1163 if (present(time_present)) then
1164 if (time_present) td_calc_mode = .true.
1165 end if
1166
1167 if (present(kick_time)) then
1168 is_time_for_kick = kick_time
1169 if (kick_time) after_kick = .true.
1170 end if
1171
1172 select case (pcm%initial_asc)
1173 case (1)
1174 input_asc_e = .true.
1175 input_asc_ext = .false.
1176 case (2)
1177 input_asc_e = .false.
1178 input_asc_ext = .true.
1179 case (3)
1180 input_asc_e = .true.
1181 input_asc_ext = .true.
1182 case default
1183 input_asc_e = .false.
1184 input_asc_ext = .false.
1185 end select
1186
1187 if (calc == pcm_nuclei) then
1188 call pcm_v_nuclei_cav(pcm%v_n, ions, pcm%tess, pcm%n_tesserae)
1189 call pcm_charges(pcm%q_n, pcm%qtot_n, pcm%v_n, pcm%matrix, pcm%n_tesserae, &
1190 pcm%q_n_nominal, pcm%epsilon_0, pcm%renorm_charges, pcm%q_tot_tol, pcm%deltaQ_n)
1191 if (pcm%calc_method == pcm_calc_poisson) then
1192 call pcm_charge_density(pcm, pcm%q_n, pcm%qtot_n, mesh, pcm%rho_n)
1193 end if
1194 call pcm_pot_rs(pcm, pcm%v_n_rs, pcm%q_n, pcm%rho_n, mesh, psolver)
1195 end if
1196
1197 if (calc == pcm_electrons) then
1198 call pcm_v_cav_li(pcm%v_e, v_h, pcm, mesh)
1199 if (td_calc_mode .and. .not. is_close(pcm%epsilon_infty, pcm%epsilon_0) .and. pcm%tdlevel /= pcm_td_eq) then
1200
1201 select case (pcm%tdlevel)
1202 case (pcm_td_eom)
1203 select case (pcm%which_eps)
1204 case (pcm_drude_model)
1205 call pcm_charges_propagation(pcm%q_e, pcm%v_e, pcm%dt, pcm%tess, input_asc_e, calc, &
1206 pcm_drude_model, pcm%namespace, this_drl = pcm%drl)
1207 case default
1208 call pcm_charges_propagation(pcm%q_e, pcm%v_e, pcm%dt, pcm%tess, input_asc_e, calc, &
1209 pcm_debye_model, pcm%namespace, this_deb = pcm%deb)
1210 end select
1211 if ((.not. pcm%localf) .and. not_yet_called) call pcm_eom_enough_initial(not_yet_called)
1212
1213 pcm%qtot_e = sum(pcm%q_e)
1214 case (pcm_td_neq)
1215 select case (pcm%iter)
1216 case (1)
1217
1220 call pcm_charges(pcm%q_e, pcm%qtot_e, pcm%v_e, pcm%matrix, pcm%n_tesserae, &
1221 pcm%q_e_nominal, pcm%epsilon_0, pcm%renorm_charges, pcm%q_tot_tol, pcm%deltaQ_e)
1222
1224 call pcm_charges(pcm%q_e_in, pcm%qtot_e_in, pcm%v_e, pcm%matrix_d, pcm%n_tesserae, &
1225 pcm%q_e_nominal, pcm%epsilon_infty, pcm%renorm_charges, pcm%q_tot_tol, pcm%deltaQ_e)
1226
1227 pcm%q_e_in = pcm%q_e - pcm%q_e_in
1228 pcm%qtot_e_in = pcm%qtot_e - pcm%qtot_e_in
1229 case default
1230
1231 call pcm_charges(pcm%q_e, pcm%qtot_e, pcm%v_e, pcm%matrix_d, pcm%n_tesserae, &
1232 pcm%q_e_nominal, pcm%epsilon_infty, pcm%renorm_charges, pcm%q_tot_tol, pcm%deltaQ_e)
1233
1234 pcm%q_e = pcm%q_e + pcm%q_e_in
1235 pcm%qtot_e = pcm%qtot_e + pcm%qtot_e_in
1236 end select
1237 end select
1238 else
1239
1240
1241 call pcm_charges(pcm%q_e, pcm%qtot_e, pcm%v_e, pcm%matrix, pcm%n_tesserae, &
1242 pcm%q_e_nominal, pcm%epsilon_0, pcm%renorm_charges, pcm%q_tot_tol, pcm%deltaQ_e)
1243
1244 end if
1245 if (pcm%calc_method == pcm_calc_poisson) then
1246 call pcm_charge_density(pcm, pcm%q_e, pcm%qtot_e, mesh, pcm%rho_e)
1247 end if
1248 call pcm_pot_rs(pcm, pcm%v_e_rs, pcm%q_e, pcm%rho_e, mesh, psolver)
1249 end if
1250
1251 if ((calc == pcm_external_plus_kick .or. calc == pcm_kick) .and. .not. pcm%kick_like) then
1252 pcm%q_ext = m_zero
1253 pcm%v_ext_rs = m_zero
1254 end if
1255
1256 if (calc == pcm_external_potential .or. calc == pcm_external_plus_kick) then
1257 call pcm_v_cav_li(pcm%v_ext, v_ext, pcm, mesh)
1258 if (td_calc_mode .and. .not. is_close(pcm%epsilon_infty, pcm%epsilon_0) .and. pcm%tdlevel /= pcm_td_eq) then
1259
1260 select case (pcm%tdlevel)
1261 case (pcm_td_eom)
1262 select case (pcm%which_eps)
1263 case (pcm_drude_model)
1264 call pcm_charges_propagation(pcm%q_ext, pcm%v_ext, pcm%dt, pcm%tess, input_asc_ext, calc, &
1265 pcm%which_eps, pcm%namespace, this_drl=pcm%drl)
1266 case default
1267 call pcm_charges_propagation(pcm%q_ext, pcm%v_ext, pcm%dt, pcm%tess, input_asc_ext, calc, &
1268 pcm%which_eps, pcm%namespace, this_deb=pcm%deb)
1269 end select
1270 if ((.not. pcm%kick_is_present) .and. not_yet_called) call pcm_eom_enough_initial(not_yet_called)
1271
1272 pcm%qtot_ext = sum(pcm%q_ext)
1273
1275 pcm%q_ext = pcm%q_ext + pcm%q_ext_in
1276 pcm%qtot_ext = pcm%qtot_ext + pcm%qtot_ext_in
1277 case (pcm_td_neq)
1279 call pcm_charges(pcm%q_ext, pcm%qtot_ext, pcm%v_ext, pcm%matrix_lf_d, pcm%n_tesserae)
1280
1281
1283 pcm%q_ext = pcm%q_ext + pcm%q_ext_in
1284 pcm%qtot_ext = pcm%qtot_ext + pcm%qtot_ext_in
1285 end select
1286 else
1287
1288
1291 call pcm_charges(pcm%q_ext, pcm%qtot_ext, pcm%v_ext, pcm%matrix_lf, pcm%n_tesserae)
1292
1293
1294 pcm%q_ext_in = pcm%q_ext
1295 pcm%qtot_ext_in = pcm%qtot_ext
1296 end if
1297 if (pcm%calc_method == pcm_calc_poisson) then
1298 call pcm_charge_density(pcm, pcm%q_ext, pcm%qtot_ext, mesh, pcm%rho_ext)
1299 end if
1300 call pcm_pot_rs(pcm, pcm%v_ext_rs, pcm%q_ext, pcm%rho_ext, mesh, psolver)
1301
1302 end if
1303
1304 if (calc == pcm_external_plus_kick .or. calc == pcm_kick) then
1305 pcm%q_kick = m_zero
1306 pcm%v_kick_rs = m_zero
1307 if (is_time_for_kick) then
1308 call pcm_v_cav_li(pcm%v_kick, kick, pcm, mesh)
1309 else
1310 pcm%v_kick = m_zero
1311 end if
1312 if (pcm%kick_like) then
1313 if (is_time_for_kick) then
1314
1315 if (.not. is_close(pcm%epsilon_infty, pcm%epsilon_0)) then
1316 call pcm_charges(pcm%q_kick, pcm%qtot_kick, pcm%v_kick, pcm%matrix_lf_d, pcm%n_tesserae)
1317 else
1318 call pcm_charges(pcm%q_kick, pcm%qtot_kick, pcm%v_kick, pcm%matrix_lf, pcm%n_tesserae)
1319 end if
1320 else
1321 pcm%q_kick = m_zero
1322 pcm%qtot_kick = m_zero
1323 if (pcm%calc_method == pcm_calc_poisson) pcm%rho_kick = m_zero
1324 pcm%v_kick_rs = m_zero
1325
1326 pop_sub(pcm_calc_pot_rs)
1327 return
1328 end if
1329 else
1330 if (after_kick) then
1331
1332 select case (pcm%which_eps)
1333 case (pcm_drude_model)
1334 call pcm_charges_propagation(pcm%q_kick, pcm%v_kick, pcm%dt, pcm%tess, input_asc_ext, calc, &
1335 pcm%which_eps, pcm%namespace, this_drl=pcm%drl)
1336 case default
1337 call pcm_charges_propagation(pcm%q_kick, pcm%v_kick, pcm%dt, pcm%tess, input_asc_ext, calc, &
1338 pcm%which_eps, pcm%namespace, this_deb=pcm%deb)
1339 end select
1340 if (not_yet_called) call pcm_eom_enough_initial(not_yet_called)
1341 pcm%qtot_kick = sum(pcm%q_kick)
1342
1343 else
1344 pcm%q_kick = m_zero
1345 pcm%qtot_kick = m_zero
1346 if (pcm%calc_method == pcm_calc_poisson) pcm%rho_kick = m_zero
1347 pcm%v_kick_rs = m_zero
1348
1349 pop_sub(pcm_calc_pot_rs)
1350 return
1351 end if
1352 end if
1353
1354 if (pcm%calc_method == pcm_calc_poisson) then
1355 call pcm_charge_density(pcm, pcm%q_kick, pcm%qtot_kick, mesh, pcm%rho_kick)
1356 end if
1357 call pcm_pot_rs(pcm, pcm%v_kick_rs, pcm%q_kick, pcm%rho_kick, mesh, psolver)
1358 if (.not. pcm%kick_like) then
1359 if (calc == pcm_external_plus_kick) then
1360 pcm%q_ext = pcm%q_ext + pcm%q_kick
1361 pcm%v_ext_rs = pcm%v_ext_rs + pcm%v_kick_rs
1362 else
1363 pcm%q_ext = pcm%q_kick
1364 pcm%v_ext_rs = pcm%v_kick_rs
1365 end if
1366 end if
1367
1368 end if
1369
1370
1371 if (mpi_world%is_root()) then
1372 write (asc_vs_t_unit_format_tail,'(I5,A11)') pcm%n_tesserae,'(1X,F14.8))'
1373 write (asc_vs_t_unit_format,'(A)') '(F14.8,'//trim(adjustl(asc_vs_t_unit_format_tail))
1374
1375 if (pcm%solute .and. pcm%localf .and. td_calc_mode .and. calc == pcm_electrons) then
1376 asc_vs_t_unit_e = io_open(pcm_dir//'ASC_e_vs_t.dat', pcm%namespace, &
1377 action='write', position='append', form='formatted')
1378 write(asc_vs_t_unit_e,trim(adjustl(asc_vs_t_unit_format))) pcm%iter*pcm%dt, &
1379 (pcm%q_e(ia) , ia = 1,pcm%n_tesserae)
1380 call io_close(asc_vs_t_unit_e)
1381 end if
1382 end if
1383
1384 pop_sub(pcm_calc_pot_rs)
1385 end subroutine pcm_calc_pot_rs
1386
1387 ! -----------------------------------------------------------------------------
1390 subroutine pcm_v_cav_li(v_cav, v_mesh, pcm, mesh)
1391 type(pcm_t), intent(in) :: pcm
1392 type(mesh_t), intent(in) :: mesh
1393 real(real64), intent(in) :: v_mesh(:)
1394 real(real64), intent(out) :: v_cav(:)
1395
1396 integer :: ia
1397
1398 type(mesh_interpolation_t) :: mesh_interpolation
1399
1400 push_sub(pcm_v_cav_li)
1401
1402 v_cav = m_zero
1403
1404 call mesh_interpolation_init(mesh_interpolation, mesh)
1405
1406 do ia = 1, pcm%n_tesserae
1407 call mesh_interpolation_evaluate(mesh_interpolation, v_mesh, pcm%tess(ia)%point, v_cav(ia))
1408 end do
1409
1410 call mesh_interpolation_end(mesh_interpolation)
1411
1412 pop_sub(pcm_v_cav_li)
1413 end subroutine pcm_v_cav_li
1414
1415 !--------------------------------------------------------------------------------
1416
1419 subroutine pcm_v_nuclei_cav(v_n_cav, ions, tess, n_tess)
1420 real(real64), intent(out) :: v_n_cav(:)
1421 type(ions_t), intent(in) :: ions
1422 type(pcm_tessera_t), intent(in) :: tess(:)
1423 integer, intent(in) :: n_tess
1424
1425 real(real64) :: dist
1426 integer :: ik, ia
1427
1428 push_sub(pcm_v_nuclei_cav)
1429
1430 v_n_cav = m_zero
1431
1432 do ik = 1, n_tess
1433 do ia = 1, ions%natoms
1434
1435 dist = norm2(ions%pos(1:pcm_dim_space, ia) - tess(ik)%point)
1436
1437 v_n_cav(ik) = v_n_cav(ik) - ions%charge(ia)/dist
1438 end do
1439 end do
1440
1441 pop_sub(pcm_v_nuclei_cav)
1442 end subroutine pcm_v_nuclei_cav
1443
1444 ! -----------------------------------------------------------------------------
1445
1449 subroutine pcm_elect_energy(ions, pcm, E_int_ee, E_int_en, E_int_ne, E_int_nn, E_int_e_ext, E_int_n_ext)
1450 type(ions_t), intent(in) :: ions
1451 type(pcm_t), intent(in) :: pcm
1452 real(real64), intent(out) :: e_int_ee
1453 real(real64), intent(out) :: e_int_en
1454 real(real64), intent(out) :: e_int_ne
1455 real(real64), intent(out) :: e_int_nn
1456 real(real64), optional, intent(out) :: e_int_e_ext
1457 real(real64), optional, intent(out) :: e_int_n_ext
1458
1459 real(real64) :: dist, z_ia
1460 integer :: ik, ia
1461
1462 push_sub(pcm_elect_energy)
1463
1464 e_int_ee = m_zero
1465 e_int_en = m_zero
1466 e_int_ne = m_zero
1467 e_int_nn = m_zero
1468
1469 if (pcm%localf .and. ((.not. present(e_int_e_ext)) .or. &
1470 (.not. present(e_int_n_ext)))) then
1471 message(1) = "pcm_elect_energy: There are lacking terms in subroutine call."
1472 call messages_fatal(1, namespace=pcm%namespace)
1473 else if (pcm%localf .and. (present(e_int_e_ext) .and. &
1474 present(e_int_n_ext))) then
1475 e_int_e_ext = m_zero
1476 e_int_n_ext = m_zero
1477 end if
1478
1479 do ik = 1, pcm%n_tesserae
1480
1481 e_int_ee = e_int_ee + pcm%v_e(ik)*pcm%q_e(ik)
1482 e_int_en = e_int_en + pcm%v_e(ik)*pcm%q_n(ik)
1483 if (pcm%localf) then
1484 e_int_e_ext = e_int_e_ext + pcm%v_e(ik)*pcm%q_ext(ik)
1485 end if
1486
1487 do ia = 1, ions%natoms
1488
1489 dist = norm2(ions%pos(1:pcm_dim_space, ia) - pcm%tess(ik)%point)
1490
1491 z_ia = -ions%charge(ia)
1492
1493 e_int_ne = e_int_ne + z_ia*pcm%q_e(ik) / dist
1494 e_int_nn = e_int_nn + z_ia*pcm%q_n(ik) / dist
1495 if (pcm%localf) e_int_n_ext = e_int_n_ext + z_ia*pcm%q_ext(ik) / dist
1496
1497 end do
1498 end do
1499
1500 e_int_ee = m_half*e_int_ee
1501 e_int_en = m_half*e_int_en
1502 e_int_ne = m_half*e_int_ne
1503 e_int_nn = m_half*e_int_nn
1504 ! E_int_e_ext and E_int_n_ext do not need 1/2 factor
1505
1506 ! print results of the iteration in pcm_info file
1507 if (mpi_world%is_root()) then
1508 if (pcm%localf) then
1509 write(pcm%info_unit, &
1510 '(3X,I5,5X,F20.8,5X,F20.8,5X,F20.8,5X,F20.8,5X,F20.8,5X,F20.8,5X,F20.8,5X,F20.8,F20.8,5X,F20.8,5X,F20.8,5X,F20.8,5X)') &
1511 pcm%counter, &
1512 units_from_atomic(units_out%energy, e_int_ee), &
1513 units_from_atomic(units_out%energy, e_int_en), &
1514 units_from_atomic(units_out%energy, e_int_nn), &
1515 units_from_atomic(units_out%energy, e_int_ne), &
1516 units_from_atomic(units_out%energy, e_int_e_ext), &
1517 units_from_atomic(units_out%energy, e_int_n_ext), &
1518 units_from_atomic(units_out%energy, e_int_ee + e_int_en + e_int_nn + e_int_ne + e_int_e_ext + e_int_n_ext), &
1519 pcm%qtot_e, &
1520 pcm%deltaQ_e, &
1521 pcm%qtot_n, &
1522 pcm%deltaQ_n, &
1523 pcm%qtot_ext
1524 else
1525 write(pcm%info_unit, &
1526 '(3X,I5,5X,F20.8,5X,F20.8,5X,F20.8,5X,F20.8,5X,F20.8,5X,F20.8,5X,F20.8,5X,F20.8,5X,F20.8)') &
1527 pcm%counter, &
1528 units_from_atomic(units_out%energy, e_int_ee), &
1529 units_from_atomic(units_out%energy, e_int_en), &
1530 units_from_atomic(units_out%energy, e_int_nn), &
1531 units_from_atomic(units_out%energy, e_int_ne), &
1532 units_from_atomic(units_out%energy, e_int_ee + e_int_en + e_int_nn + e_int_ne), &
1533 pcm%qtot_e, &
1534 pcm%deltaQ_e, &
1535 pcm%qtot_n, &
1536 pcm%deltaQ_n
1537 end if
1538 end if
1539
1540 pop_sub(pcm_elect_energy)
1541 end subroutine pcm_elect_energy
1542
1543 ! -----------------------------------------------------------------------------
1544
1548 subroutine pcm_charges(q_pcm, q_pcm_tot, v_cav, pcm_mat, n_tess, qtot_nominal, epsilon, renorm_charges, q_tot_tol, deltaQ)
1549 real(real64), intent(out) :: q_pcm(:)
1550 real(real64), intent(out) :: q_pcm_tot
1551 real(real64), intent(in) :: v_cav(:)
1552 real(real64), intent(in) :: pcm_mat(:,:)
1553 integer, intent(in) :: n_tess
1554 real(real64), optional, intent(in) :: qtot_nominal
1555 real(real64), optional, intent(in) :: epsilon
1556 logical, optional, intent(in) :: renorm_charges
1557 real(real64), optional, intent(in) :: q_tot_tol
1558 real(real64), optional, intent(out) :: deltaq
1559
1560 integer :: ia
1561 real(real64) :: q_pcm_tot_norm
1562 real(real64) :: coeff
1563
1564 push_sub(pcm_charges)
1565 call profiling_in('PCM_CHARGES')
1566
1567 ! q_pcm = pcm_mat * v_cav (beta = 0 overwrites q_pcm)
1568 call lalg_gemv(n_tess, n_tess, m_one, pcm_mat, v_cav, m_zero, q_pcm)
1569 q_pcm_tot = sum(q_pcm(1:n_tess))
1570
1571 if (present(qtot_nominal) .and. present(epsilon) .and. &
1572 present(renorm_charges) .and. present(q_tot_tol) .and. present(deltaq)) then
1573
1574 deltaq = abs(q_pcm_tot) - ((epsilon - m_one)/epsilon) * abs(qtot_nominal)
1575 if ((renorm_charges) .and. (abs(deltaq) > q_tot_tol)) then
1576 q_pcm_tot_norm = m_zero
1577 coeff = sign(m_one, qtot_nominal)*sign(m_one, deltaq)
1578 do ia = 1, n_tess
1579 q_pcm(ia) = q_pcm(ia) + coeff*q_pcm(ia)/q_pcm_tot*abs(deltaq)
1580 q_pcm_tot_norm = q_pcm_tot_norm + q_pcm(ia)
1581 end do
1582 q_pcm_tot = q_pcm_tot_norm
1583 end if
1584 end if
1585
1586 call profiling_out('PCM_CHARGES')
1587 pop_sub(pcm_charges)
1588 end subroutine pcm_charges
1589
1590 ! -----------------------------------------------------------------------------
1592 logical function pcm_nn_in_mesh(pcm, mesh) result(in_mesh)
1593 type(pcm_t), intent(in) :: pcm
1594 class(mesh_t), intent(in) :: mesh
1595
1596 integer :: ia, nm(pcm%space%dim), ipt, i1, i2, i3
1597 real(real64) :: posrel(pcm%space%dim)
1598 integer(int64) :: pt
1599
1600 push_sub(pcm_nn_in_mesh)
1601
1602 in_mesh = .true.
1603 do ia = 1, pcm%n_tesserae
1604
1605 posrel(:) = pcm%tess(ia)%point(1:pcm%space%dim)/mesh%spacing(1:pcm%space%dim)
1606
1607 nm(:) = floor(posrel(:))
1608
1609 ! Get the nearest neighboring points
1610 ipt = 0
1611 do i1 = -pcm%tess_nn + 1 , pcm%tess_nn
1612 do i2 = -pcm%tess_nn + 1 , pcm%tess_nn
1613 do i3 = -pcm%tess_nn + 1 , pcm%tess_nn
1614 ipt = ipt + 1
1615 pt = mesh_global_index_from_coords(mesh, [i1 + nm(1), i2 + nm(2), i3 + nm(3)])
1616
1617 if (pt <= 0 .or. pt > mesh%np_part_global) then
1618 in_mesh = .false.
1619 pop_sub(pcm_nn_in_mesh)
1620 return
1621 end if
1622
1623 end do
1624 end do
1625 end do
1626 end do
1627
1628 pop_sub(pcm_nn_in_mesh)
1629 end function pcm_nn_in_mesh
1630
1631 ! -----------------------------------------------------------------------------
1633 subroutine pcm_poisson_sanity_check(pcm, mesh)
1634 type(pcm_t), intent(in) :: pcm
1635 class(mesh_t), intent(in) :: mesh
1636
1637 push_sub(pcm_poisson_sanity_check)
1638
1639 if (.not. pcm_nn_in_mesh(pcm, mesh)) then
1640 message(1) = 'The simulation box is too small to contain all the requested'
1641 message(2) = 'nearest neighbors for each tessera.'
1642 message(3) = 'Consider using a larger box or reduce PCMChargeSmearNN.'
1643 call messages_warning(3, namespace=pcm%namespace)
1644 end if
1645
1646 pop_sub(pcm_poisson_sanity_check)
1647 end subroutine pcm_poisson_sanity_check
1648
1649 ! -----------------------------------------------------------------------------
1652 subroutine pcm_charge_density(pcm, q_pcm, q_pcm_tot, mesh, rho)
1653 type(pcm_t), intent(inout) :: pcm
1654 real(real64), intent(in) :: q_pcm(:)
1655 real(real64), intent(in) :: q_pcm_tot
1656 type(mesh_t), intent(in) :: mesh
1657 real(real64), intent(out) :: rho(:)
1658
1659 integer :: ia
1660 real(real64) :: norm, qtot, rr, xx(pcm%space%dim), pp(pcm%space%dim)
1661
1662 ! nearest neighbor variables
1663 integer :: nm(pcm%space%dim)
1664 real(real64) :: posrel(pcm%space%dim)
1665 integer :: npt, ipt
1666 integer :: i1, i2, i3
1667 integer, allocatable :: pt(:)
1668 real(real64), allocatable :: lrho(:) ! local charge density on a tessera NN
1669 logical :: inner_point, boundary_point
1670
1671 !profiling and debug
1672 integer :: ierr
1673
1674 push_sub(pcm_charge_density)
1675 call profiling_in('PCM_CHARGE_DENSITY')
1676
1677 npt = (2*pcm%tess_nn)**pcm%space%dim
1678 safe_allocate(pt(1:npt))
1679 safe_allocate(lrho(1:npt))
1680
1681 pt = 0
1682 rho = m_zero
1683
1684 do ia = 1, pcm%n_tesserae
1685
1686 pp(:) = pcm%tess(ia)%point(1:pcm%space%dim)
1687 posrel(:) = pp(:)/mesh%spacing(1:pcm%space%dim)
1688
1689 nm(:) = floor(posrel(:))
1690
1691 ! Get the nearest neighboring points
1692 ipt = 0
1693 do i1 = -pcm%tess_nn + 1 , pcm%tess_nn
1694 do i2 = -pcm%tess_nn + 1 , pcm%tess_nn
1695 do i3 = -pcm%tess_nn + 1 , pcm%tess_nn
1696 ipt = ipt + 1
1697 pt(ipt) = mesh_local_index_from_coords(mesh, [i1 + nm(1), i2 + nm(2), i3 + nm(3)])
1698 end do
1699 end do
1700 end do
1701
1702
1703 ! Extrapolate the tessera point charge with a gaussian distritibution
1704 ! to the neighboring points
1705 ! \f$ rho(r) = N exp(-|r-sk|^2/(alpha*Ak)) \f$
1706 norm = m_zero
1707 lrho = m_zero
1708 do ipt = 1, npt
1709
1710 ! Check the point is inside the mesh skip otherwise
1711 if (pt(ipt) > 0 .and. pt(ipt) <= mesh%np_part) then
1712
1713 if (mesh%parallel_in_domains) then
1714 boundary_point = pt(ipt) > mesh%np + mesh%pv%np_ghost
1715 inner_point = pt(ipt) > 0 .and. pt(ipt) <= mesh%np
1716
1717 if (boundary_point .or. inner_point) then
1718 xx(:) = mesh%x(:, pt(ipt))
1719 else
1720 cycle
1721 end if
1722
1723 else
1724 xx(:) = mesh%x(:, pt(ipt))
1725 end if
1726
1727 rr = sum((xx(1:pcm%space%dim) - pp(1:pcm%space%dim))**2)
1728 norm = norm + exp(-rr/(pcm%tess(ia)%area*pcm%gaussian_width))
1729 lrho(ipt) = lrho(ipt) + exp(-rr/(pcm%tess(ia)%area*pcm%gaussian_width))
1730
1731 end if
1732 end do
1733
1734 ! reduce the local density scattered among nodes
1735 call mesh%allreduce(lrho, npt)
1736
1737 ! normalize the integral to the tessera point charge q_pcm(ia)
1738 norm = sum(lrho(1:npt)) * mesh%volume_element
1739 if (norm > m_epsilon) then
1740 norm = q_pcm(ia)/norm
1741 else
1742 norm = m_zero
1743 end if
1744 lrho(:) = lrho(:) * norm
1745
1746 ! Add up the local density to the full charge density
1747 do ipt = 1, npt
1748
1749 if (pt(ipt) > 0 .and. pt(ipt) <= mesh%np_part_global) then
1750
1751 if (mesh%parallel_in_domains) then
1752 boundary_point = pt(ipt) > mesh%np + mesh%pv%np_ghost
1753 inner_point = pt(ipt) > 0 .and. pt(ipt) <= mesh%np
1754 if (boundary_point .or. inner_point) rho(pt(ipt)) = rho(pt(ipt)) + lrho(ipt)
1755 else
1756 rho(pt(ipt)) = rho(pt(ipt)) + lrho(ipt)
1757 end if
1758
1759 end if
1760 end do
1761
1762 end do
1763
1764 if (debug%info) then
1765 qtot = dmf_integrate(mesh, rho)
1766 call messages_write(' PCM charge density integrates to q = ')
1767 call messages_write(qtot)
1768 call messages_write(' (qtot = ')
1769 call messages_write(q_pcm_tot)
1770 call messages_write(')')
1771 call messages_info()
1772
1773 ! Keep this here for debug purposes.
1774 call dio_function_output(io_function_fill_how("VTK"), ".", "rho_pcm", pcm%namespace, pcm%space, &
1775 mesh, rho, unit_one, ierr)
1776 end if
1777
1778 safe_deallocate_a(pt)
1779 safe_deallocate_a(lrho)
1780
1781 call profiling_out('PCM_CHARGE_DENSITY')
1782
1783 pop_sub(pcm_charge_density)
1784 end subroutine pcm_charge_density
1785
1786
1787 ! -----------------------------------------------------------------------------
1789 subroutine pcm_pot_rs(pcm, v_pcm, q_pcm, rho, mesh, psolver)
1790 type(pcm_t), intent(inout) :: pcm
1791 real(real64), contiguous, intent(inout) :: v_pcm(:)
1792 real(real64), contiguous, intent(in) :: q_pcm(:)
1793 real(real64), contiguous, intent(inout) :: rho(:)
1794 type(mesh_t), intent(in) :: mesh
1795 type(poisson_t), intent(in) :: psolver
1796
1797 integer :: ierr
1798
1799 push_sub(pcm_pot_rs)
1800 call profiling_in('PCM_POT_RS')
1801
1802 v_pcm = m_zero
1803
1804 select case (pcm%calc_method)
1805 case (pcm_calc_direct)
1806 call pcm_pot_rs_direct(v_pcm, q_pcm, pcm%tess, pcm%n_tesserae, mesh, pcm%gaussian_width)
1807
1808 case (pcm_calc_poisson)
1809 call pcm_pot_rs_poisson(pcm%namespace, v_pcm, psolver, rho)
1810
1811 case default
1812 message(1) = "BAD BAD BAD"
1813 call messages_fatal(1,only_root_writes = .true., namespace=pcm%namespace)
1814
1815 end select
1816
1817 if (debug%info) then
1818 ! Keep this here for debug purposes.
1819 call dio_function_output(io_function_fill_how("VTK"), ".", "v_pcm", pcm%namespace, pcm%space, &
1820 mesh, v_pcm, unit_one, ierr)
1821 end if
1822
1823 call profiling_out('PCM_POT_RS')
1824 pop_sub(pcm_pot_rs)
1825 end subroutine pcm_pot_rs
1826
1827
1828 ! -----------------------------------------------------------------------------
1830 subroutine pcm_pot_rs_poisson(namespace, v_pcm, psolver, rho)
1831 type(namespace_t), intent(in) :: namespace
1832 real(real64), contiguous, intent(inout) :: v_pcm(:)
1833 type(poisson_t), intent(in) :: psolver
1834 real(real64), contiguous, intent(inout) :: rho(:)
1835
1836 push_sub(pcm_pot_rs_poisson)
1837
1838 call dpoisson_solve(psolver, namespace, v_pcm, rho)
1839
1840 pop_sub(pcm_pot_rs_poisson)
1841 end subroutine pcm_pot_rs_poisson
1842
1843
1844 ! -----------------------------------------------------------------------------
1846 subroutine pcm_pot_rs_direct(v_pcm, q_pcm, tess, n_tess, mesh, width_factor)
1847 real(real64), intent(out) :: v_pcm(:)
1848 real(real64), intent(in) :: q_pcm(:)
1849 real(real64), intent(in) :: width_factor
1850 integer, intent(in) :: n_tess
1851 type(mesh_t), intent(in) :: mesh
1852 type(pcm_tessera_t), intent(in) :: tess(:)
1853
1854 real(real64), parameter :: p_1 = 0.119763_real64
1855 real(real64), parameter :: p_2 = 0.205117_real64
1856 real(real64), parameter :: q_1 = 0.137546_real64
1857 real(real64), parameter :: q_2 = 0.434344_real64
1858 real(real64) :: arg
1859 real(real64) :: term
1860 integer :: ip
1861 integer :: ia
1862
1863 push_sub(pcm_pot_rs_direct)
1864
1865 v_pcm = m_zero
1866
1867 if (abs(width_factor) > m_epsilon) then
1868
1869 do ia = 1, n_tess
1870 do ip = 1, mesh%np
1871 ! Computing the distances between tesserae and grid points.
1872 call mesh_r(mesh, ip, term, origin=tess(ia)%point)
1873 arg = term/sqrt(tess(ia)%area*width_factor)
1874 term = (m_one + p_1*arg + p_2*arg**2) / (m_one + q_1*arg + q_2*arg**2 + p_2*arg**3)
1875 v_pcm(ip) = v_pcm(ip) + q_pcm(ia) * term/sqrt(tess(ia)%area*width_factor)
1876 end do
1877 end do
1878
1879 v_pcm = m_two*v_pcm / sqrt(m_pi)
1880
1881 else
1882
1883 do ia = 1, n_tess
1884 do ip = 1, mesh%np
1885 ! Computing the distances between tesserae and grid points.
1886 call mesh_r(mesh, ip, term, origin=tess(ia)%point)
1887 v_pcm(ip) = v_pcm(ip) + q_pcm(ia)/term
1888 end do
1889 end do
1890 end if
1891
1892 pop_sub(pcm_pot_rs_direct)
1893 end subroutine pcm_pot_rs_direct
1894
1895 ! -----------------------------------------------------------------------------
1896
1898 subroutine pcm_matrix(eps, tess, n_tess, pcm_mat, localf)
1899 real(real64), intent(in) :: eps
1900 type(pcm_tessera_t), intent(in) :: tess(:)
1901 integer, intent(in) :: n_tess
1902 real(real64), intent(out) :: pcm_mat(:,:)
1903 logical, optional, intent(in) :: localf
1904
1905 integer :: i, info
1906 integer, allocatable :: iwork(:)
1907 real(real64), allocatable :: mat_tmp(:,:)
1908
1909 real(real64) :: sgn_lf
1910
1911 push_sub(pcm_matrix)
1912
1914 safe_allocate(s_mat_act(1:n_tess, 1:n_tess))
1915 call s_i_matrix(n_tess, tess)
1916
1918 safe_allocate(sigma(1:n_tess, 1:n_tess))
1919 sigma = s_mat_act/eps
1920
1922 safe_allocate(d_mat_act(1:n_tess, 1:n_tess))
1923 call d_i_matrix(n_tess, tess)
1924
1926 safe_allocate(delta(1:n_tess, 1:n_tess))
1927 delta = d_mat_act
1928
1929 sgn_lf = m_one
1931 if (present(localf)) then
1932 if (localf) sgn_lf = -m_one
1933 end if
1934
1936 pcm_mat = - sgn_lf * d_mat_act
1937
1938 do i = 1, n_tess
1939 pcm_mat(i,i) = pcm_mat(i,i) + m_two*m_pi
1940 end do
1941
1942 safe_deallocate_a(d_mat_act)
1943
1944 safe_allocate(iwork(1:n_tess))
1945
1948 ! FIXME: use interface, or routine in lalg_adv_lapack_inc.F90
1949 call dgesv(n_tess, n_tess, s_mat_act, n_tess, iwork, pcm_mat, n_tess, info)
1950
1951 safe_deallocate_a(iwork)
1952
1953 safe_deallocate_a(s_mat_act)
1954
1957 pcm_mat = -matmul(sigma, pcm_mat)
1958
1959 do i = 1, n_tess
1960 pcm_mat(i,i) = pcm_mat(i,i) + m_two*m_pi
1961 end do
1962
1963 pcm_mat = pcm_mat - sgn_lf * delta
1964
1965 safe_allocate(mat_tmp(1:n_tess,1:n_tess))
1966 mat_tmp = m_zero
1967
1968 safe_allocate(d_mat_act(1:n_tess,1:n_tess))
1969 call d_i_matrix(n_tess, tess)
1970
1971 mat_tmp = transpose(d_mat_act)
1972
1973 mat_tmp = matmul(sigma, mat_tmp)
1974
1975 mat_tmp = mat_tmp + m_two*m_pi*sigma
1976
1977 safe_deallocate_a(d_mat_act)
1978
1979 safe_allocate(s_mat_act(1:n_tess, 1:n_tess))
1980 call s_i_matrix(n_tess, tess)
1981
1982 mat_tmp = mat_tmp + m_two*m_pi*s_mat_act - matmul(delta, s_mat_act)
1983
1984 safe_deallocate_a(s_mat_act)
1985 safe_deallocate_a(sigma)
1986 safe_deallocate_a(delta)
1987
1988 safe_allocate(iwork(1:n_tess))
1989
1992 call dgesv(n_tess, n_tess, mat_tmp, n_tess, iwork, pcm_mat, n_tess, info)
1993
1994 safe_deallocate_a(iwork)
1995 safe_deallocate_a(mat_tmp)
1996
1997 pcm_mat = - sgn_lf * pcm_mat
1998
1999 ! Testing
2000 if (gamess_benchmark) then
2001 do i = 1, n_tess
2002 mat_gamess(i,:) = pcm_mat(i,:)/tess(i)%area
2003 end do
2004 end if
2005
2006 pop_sub(pcm_matrix)
2007 end subroutine pcm_matrix
2008
2009 ! -----------------------------------------------------------------------------
2010
2011 subroutine s_i_matrix(n_tess, tess)
2012 integer, intent(in) :: n_tess
2013 type(pcm_tessera_t), intent(in) :: tess(:)
2014
2015 integer :: ii, jj
2016
2017 s_mat_act = m_zero
2018
2019 do ii = 1, n_tess
2020 do jj = ii, n_tess
2021 s_mat_act(ii,jj) = s_mat_elem_i(tess(ii), tess(jj))
2022 if (ii /= jj) s_mat_act(jj,ii) = s_mat_act(ii,jj) !symmetric matrix
2023 end do
2024 end do
2025
2026 end subroutine s_i_matrix
2027
2028 ! -----------------------------------------------------------------------------
2029
2030 subroutine d_i_matrix(n_tess, tess)
2031 integer, intent(in) :: n_tess
2032 type(pcm_tessera_t), intent(in) :: tess(:)
2033
2034 integer :: ii, jj
2035
2036 d_mat_act = m_zero
2037
2038 do ii = 1, n_tess
2039 do jj = 1, n_tess
2040 d_mat_act(ii,jj) = d_mat_elem_i(tess(ii), tess(jj))
2041 end do
2042 end do
2043
2044 end subroutine d_i_matrix
2045
2046 ! -----------------------------------------------------------------------------
2047
2050 real(real64) function s_mat_elem_I(tessi, tessj)
2051 type(pcm_tessera_t), intent(in) :: tessi
2052 type(pcm_tessera_t), intent(in) :: tessj
2053
2054 real(real64), parameter :: M_SD_DIAG = 1.0694_real64
2055 real(real64), parameter :: M_DIST_MIN = 0.1_real64
2056
2057 real(real64) :: diff(1:PCM_DIM_SPACE)
2058 real(real64) :: dist
2059 real(real64) :: s_diag
2060 real(real64) :: s_off_diag
2061
2062 s_diag = m_zero
2063 s_off_diag = m_zero
2064
2065 diff = tessi%point - tessj%point
2066
2067 dist = norm2(diff)
2068
2069 if (abs(dist) <= m_epsilon) then
2070 s_diag = m_sd_diag*sqrt(m_four*m_pi/tessi%area)
2071 s_mat_elem_i = s_diag
2072 else
2073 if (dist > m_dist_min) s_off_diag = m_one/dist
2074 s_mat_elem_i = s_off_diag
2075 end if
2076
2077 end function s_mat_elem_i
2078
2079 ! -----------------------------------------------------------------------------
2080
2082 real(real64) function d_mat_elem_I(tessi, tessj)
2083 type(pcm_tessera_t), intent(in) :: tessi
2084 type(pcm_tessera_t), intent(in) :: tessj
2085
2086 real(real64), parameter :: M_SD_DIAG = 1.0694_real64
2087 real(real64), parameter :: M_DIST_MIN = 0.04_real64
2088
2089 real(real64) :: diff(1:PCM_DIM_SPACE)
2090 real(real64) :: dist
2091 real(real64) :: d_diag
2092 real(real64) :: d_off_diag
2093
2094 d_diag = m_zero
2095 d_off_diag = m_zero
2096 diff = m_zero
2097
2098 diff = tessi%point - tessj%point
2099
2100 dist = norm2(diff)
2101
2102 if (abs(dist) <= m_epsilon) then
2104 d_diag = m_sd_diag*sqrt(m_four*m_pi*tessi%area)
2105 d_diag = -d_diag/(m_two*tessi%r_sphere)
2106 d_mat_elem_i = d_diag
2107
2108 else
2110 if (dist > m_dist_min) then
2111 d_off_diag = dot_product( diff, tessj%normal(:))
2112 d_off_diag = d_off_diag*tessj%area/dist**3
2113 end if
2114
2115 d_mat_elem_i = d_off_diag
2116
2117 end if
2118
2119 end function d_mat_elem_i
2120
2121 ! -----------------------------------------------------------------------------
2122
2126 subroutine cav_gen(tess_sphere, tess_min_distance, nesf, sfe, nts, cts, unit_pcminfo)
2127 integer, intent(in) :: tess_sphere
2128 real(real64) , intent(in) :: tess_min_distance
2129 type(pcm_sphere_t), intent(inout) :: sfe(:)
2130 integer, intent(in) :: nesf
2131 integer, intent(out) :: nts
2132 type(pcm_tessera_t), intent(out) :: cts(:)
2133 integer, intent(in) :: unit_pcminfo
2134
2135 integer, parameter :: DIM_ANGLES = 24
2136 integer, parameter :: DIM_TEN = 10
2137 integer, parameter :: DIM_VERTICES = 122
2138 integer, parameter :: MAX_VERTICES = 6
2139 integer, parameter :: MXTS = 10000
2140
2141 real(real64), save :: thev(1:DIM_ANGLES)
2142 real(real64), save :: fiv(1:DIM_ANGLES)
2143 real(real64), save :: fir
2144 real(real64) :: cv(1:DIM_VERTICES, 1:PCM_DIM_SPACE)
2145 real(real64) :: th
2146 real(real64) :: fi
2147 real(real64) :: cth
2148 real(real64) :: sth
2149
2150 real(real64) :: xctst(tess_sphere*n_tess_sphere)
2151 real(real64) :: yctst(tess_sphere*n_tess_sphere)
2152 real(real64) :: zctst(tess_sphere*n_tess_sphere)
2153 real(real64) :: ast(tess_sphere*n_tess_sphere)
2154 real(real64) :: nctst(pcm_dim_space, tess_sphere*n_tess_sphere)
2155
2156 real(real64) :: pts(1:pcm_dim_space, 1:dim_ten)
2157 real(real64) :: pp(1:pcm_dim_space)
2158 real(real64) :: pp1(1:pcm_dim_space)
2159 real(real64) :: ccc(1:pcm_dim_space, 1:dim_ten)
2160
2161 integer, save :: idum(1:n_tess_sphere*max_vertices)
2162 integer :: jvt1(1:max_vertices,1:n_tess_sphere)
2163 integer :: isfet(1:dim_ten*dim_angles)
2164
2165 integer :: ii
2166 integer :: ia
2167 integer :: ja
2168 integer :: nn
2169 integer :: nsfe
2170 integer :: its
2171 integer :: n1
2172 integer :: n2
2173 integer :: n3
2174 integer :: nv
2175 integer :: i_tes
2176
2177 real(real64) :: xen
2178 real(real64) :: yen
2179 real(real64) :: zen
2180 real(real64) :: ren
2181 real(real64) :: area
2182 real(real64) :: test2
2183 real(real64) :: rij
2184 real(real64) :: dnorm
2185
2186 real(real64) :: xi
2187 real(real64) :: yi
2188 real(real64) :: zi
2189 real(real64) :: xj
2190 real(real64) :: yj
2191 real(real64) :: zj
2192
2193 real(real64) :: vol
2194 real(real64) :: stot
2195 real(real64) :: prod
2196 real(real64) :: dr
2197
2198 logical :: band_iter
2199
2200 push_sub(cav_gen)
2201
2204 data thev/ 0.6523581398_real64 , 1.107148718_real64 , 1.382085796_real64 , &
2205 1.759506858_real64 , 2.034443936_real64 , 2.489234514_real64 , &
2206 0.3261790699_real64 , 0.5535743589_real64, &
2207 0.8559571251_real64 , 0.8559571251_real64 , 1.017221968_real64 , &
2208 1.229116717_real64 , 1.229116717_real64 , 1.433327788_real64 , &
2209 1.570796327_real64 , 1.570796327_real64 , 1.708264866_real64 , &
2210 1.912475937_real64 , 1.912475937_real64 , 2.124370686_real64 , &
2211 2.285635528_real64 , 2.285635528_real64 , 2.588018295_real64 , &
2212 2.815413584_real64 /
2213 data fiv/ 0.6283185307_real64 , m_zero , &
2214 0.6283185307_real64 , m_zero , 0.6283185307_real64, &
2215 m_zero , 0.6283185307_real64 , m_zero, &
2216 0.2520539002_real64 , 1.004583161_real64 , 0.6283185307_real64, &
2217 0.3293628477_real64 , 0.9272742138_real64 , m_zero, &
2218 0.3141592654_real64 , 0.9424777961_real64 , 0.6283185307_real64, &
2219 0.2989556830_real64 , 0.9576813784_real64 , m_zero, &
2220 0.3762646305_real64 , 0.8803724309_real64 , 0.6283188307_real64, &
2221 m_zero /
2222 data fir / 1.256637061_real64 /
2223
2226 data (idum(ii),ii = 1, 280) / &
2227 1, 6, 2, 32, 36, 37, 1, 2, 3, 33, 32, 38, 1, 3, 4, 34, &
2228 33, 39, 1, 4, 5, 35, 34, 40, 1, 5, 6, 36, 35, 41, 7, 2, 6, 51, &
2229 42, 37, 8, 3, 2, 47, 43, 38, 9, 4, 3, 48, 44, 39, 10, 5, 4, &
2230 49, 45, 40, 11, 6, 5, 50, 46, 41, 8, 2, 12, 62, 47, 52, 9, &
2231 3, 13, 63, 48, 53, 10, 4, 14, 64, 49, 54, 11, 5, 15, 65, 50, &
2232 55, 7, 6, 16, 66, 51, 56, 7, 12, 2, 42, 57, 52, 8, 13, 3, &
2233 43, 58, 53, 9, 14, 4, 44, 59, 54, 10, 15, 5, 45, 60, 55, 11, &
2234 16, 6, 46, 61, 56, 8, 12, 18, 68, 62, 77, 9, 13, 19, 69, 63, &
2235 78, 10, 14, 20, 70, 64, 79, 11, 15, 21, 71, 65, 80, 7, 16, &
2236 17, 67, 66, 81, 7, 17, 12, 57, 67, 72, 8, 18, 13, 58, 68, 73, &
2237 9, 19, 14, 59, 69, 74, 10, 20, 15, 60, 70, 75, 11, 21, 16, &
2238 61, 71, 76, 22, 12, 17, 87, 82, 72, 23, 13, 18, 88, 83, 73, &
2239 24, 14, 19, 89, 84, 74, 25, 15, 20, 90, 85, 75, 26, 16, 21, &
2240 91, 86, 76, 22, 18, 12, 82, 92, 77, 23, 19, 13, 83, 93, 78, &
2241 24, 20, 14, 84, 94, 79, 25, 21, 15, 85, 95, 80, 26, 17, 16, &
2242 86, 96, 81, 22, 17, 27, 102, 87, 97, 23, 18, 28, 103, 88, 98, &
2243 24, 19, 29, 104, 89, 99, 25, 20, 30, 105, 90, 100, 26, 21, &
2244 31, 106, 91, 101, 22, 28, 18, 92, 107, 98, 23, 29, 19, 93 /
2245 data (idum(ii),ii = 281,360) / &
2246 108, 99, 24, 30, 20, 94, 109, 100, 25, 31, 21, 95, 110, 101, &
2247 26, 27, 17, 96, 111, 97, 22, 27, 28, 107, 102, 112, 23, 28, &
2248 29, 108, 103, 113, 24, 29, 30, 109, 104, 114, 25, 30, 31, &
2249 110, 105, 115, 26, 31, 27, 111, 106, 116, 122, 28, 27, 117, &
2250 118, 112, 122, 29, 28, 118, 119, 113, 122, 30, 29, 119, 120, &
2251 114, 122, 31, 30, 120, 121, 115, 122, 27, 31, 121, 117, 116 /
2252
2253 if (mpi_world%is_root()) then
2254 if (tess_sphere == 1) then
2255 write(unit_pcminfo,'(A1)') '#'
2256 write(unit_pcminfo,'(A34)') '# Number of tesserae / sphere = 60'
2257 write(unit_pcminfo,'(A1)') '#'
2258 else
2259 write(unit_pcminfo,'(A1)') '#'
2260 write(unit_pcminfo,'(A35)') '# Number of tesserae / sphere = 240'
2261 write(unit_pcminfo,'(A1)') '#'
2262 end if
2263 end if
2264
2267 dr = 0.01_real64
2268 dr = dr*p_a_b
2269
2270 sfe(:)%x = sfe(:)%x*p_a_b
2271 sfe(:)%y = sfe(:)%y*p_a_b
2272 sfe(:)%z = sfe(:)%z*p_a_b
2273 sfe(:)%r = sfe(:)%r*p_a_b
2274
2275 vol = m_zero
2276 stot = m_zero
2277 jvt1 = reshape(idum,(/6,60/))
2278
2290
2291 cv(1,1) = m_zero
2292 cv(1,2) = m_zero
2293 cv(1,3) = m_one
2294
2295 cv(122,1) = m_zero
2296 cv(122,2) = m_zero
2297 cv(122,3) = -m_one
2298
2299 ii = 1
2300 do ia = 1, dim_angles
2301 th = thev(ia)
2302 fi = fiv(ia)
2303 cth = cos(th)
2304 sth = sin(th)
2305 do ja = 1, 5
2306 fi = fi + fir
2307 if (ja == 1) fi = fiv(ia)
2308 ii = ii + 1
2309 cv(ii,1) = sth*cos(fi)
2310 cv(ii,2) = sth*sin(fi)
2311 cv(ii,3) = cth
2312 end do
2313 end do
2314
2316 nn = 0
2317 do nsfe = 1, nesf
2318 xen = sfe(nsfe)%x
2319 yen = sfe(nsfe)%y
2320 zen = sfe(nsfe)%z
2321 ren = sfe(nsfe)%r
2322
2323 xctst(:) = m_zero
2324 yctst(:) = m_zero
2325 zctst(:) = m_zero
2326 ast(:) = m_zero
2327
2328 do its = 1, n_tess_sphere
2329 do i_tes = 1, tess_sphere
2330 if (tess_sphere == 1) then
2331 n1 = jvt1(1,its)
2332 n2 = jvt1(2,its)
2333 n3 = jvt1(3,its)
2334 else
2335 if (i_tes == 1) then
2336 n1 = jvt1(1,its)
2337 n2 = jvt1(5,its)
2338 n3 = jvt1(4,its)
2339 elseif (i_tes == 2) then
2340 n1 = jvt1(4,its)
2341 n2 = jvt1(6,its)
2342 n3 = jvt1(3,its)
2343 elseif (i_tes == 3) then
2344 n1 = jvt1(4,its)
2345 n2 = jvt1(5,its)
2346 n3 = jvt1(6,its)
2347 elseif (i_tes == 4) then
2348 n1 = jvt1(2,its)
2349 n2 = jvt1(6,its)
2350 n3 = jvt1(5,its)
2351 end if
2352 end if
2353
2354 pts(1,1) = cv(n1,1)*ren + xen
2355 pts(2,1) = cv(n1,3)*ren + yen
2356 pts(3,1) = cv(n1,2)*ren + zen
2357
2358 pts(1,2) = cv(n2,1)*ren + xen
2359 pts(2,2) = cv(n2,3)*ren + yen
2360 pts(3,2) = cv(n2,2)*ren + zen
2361
2362 pts(1,3) = cv(n3,1)*ren + xen
2363 pts(2,3) = cv(n3,3)*ren + yen
2364 pts(3,3) = cv(n3,2)*ren + zen
2365
2366 pp(:) = m_zero
2367 pp1(:) = m_zero
2368 nv = 3
2369
2370 call subtessera(sfe, nsfe, nesf, nv, pts ,ccc, pp, pp1, area)
2371
2372 if (abs(area) <= m_epsilon) cycle
2373
2374 xctst(tess_sphere*(its-1) + i_tes) = pp(1)
2375 yctst(tess_sphere*(its-1) + i_tes) = pp(2)
2376 zctst(tess_sphere*(its-1) + i_tes) = pp(3)
2377 nctst(:,tess_sphere*(its-1) + i_tes) = pp1(:)
2378 ast(tess_sphere*(its-1) + i_tes) = area
2379 isfet(tess_sphere*(its-1) + i_tes) = nsfe
2380
2381 end do
2382 end do
2383
2384 do its = 1, n_tess_sphere*tess_sphere
2385
2386 if (abs(ast(its)) <= m_epsilon) cycle
2387 nn = nn + 1
2388
2389 if (nn > mxts) then
2390 write(message(1),'(a,I5,a,I5)') "total number of tesserae", nn, ">",mxts
2391 call messages_warning(1)
2392 end if
2393
2394 cts(nn)%point(1) = xctst(its)
2395 cts(nn)%point(2) = yctst(its)
2396 cts(nn)%point(3) = zctst(its)
2397 cts(nn)%normal(:) = nctst(:,its)
2398 cts(nn)%area = ast(its)
2399 cts(nn)%r_sphere = sfe(isfet(its))%r
2400
2401 end do
2402 end do
2403
2404 nts = nn
2405
2407 test2 = tess_min_distance*tess_min_distance
2408
2409 band_iter = .false.
2410 do while (.not.(band_iter))
2411 band_iter = .true.
2412
2413 loop_ia: do ia = 1, nts-1
2414 if (abs(cts(ia)%area) <= m_epsilon) cycle
2415 xi = cts(ia)%point(1)
2416 yi = cts(ia)%point(2)
2417 zi = cts(ia)%point(3)
2418
2419 loop_ja: do ja = ia+1, nts
2420 if (abs(cts(ja)%area) <= m_epsilon) cycle
2421 xj = cts(ja)%point(1)
2422 yj = cts(ja)%point(2)
2423 zj = cts(ja)%point(3)
2424
2425 rij = (xi - xj)**2 + (yi - yj)**2 + (zi - zj)**2
2426
2427 if (rij > test2) cycle
2428
2429 if (mpi_world%is_root()) then
2430 write(unit_pcminfo,'(A40,I4,A5,I4,A4,F8.4,A13,F8.4,A3)') &
2431 '# Warning: The distance between tesserae', &
2432 ia,' and ', ja,' is ',sqrt(rij),' A, less than', tess_min_distance,' A.'
2433 end if
2434
2436 xi = (xi*cts(ia)%area + xj*cts(ja)%area) / (cts(ia)%area + cts(ja)%area)
2437 yi = (yi*cts(ia)%area + yj*cts(ja)%area) / (cts(ia)%area + cts(ja)%area)
2438 zi = (zi*cts(ia)%area + zj*cts(ja)%area) / (cts(ia)%area + cts(ja)%area)
2439
2440 cts(ia)%point(1) = xi
2441 cts(ia)%point(2) = yi
2442 cts(ia)%point(3) = zi
2443
2445 cts(ia)%normal = (cts(ia)%normal*cts(ia)%area + cts(ja)%normal*cts(ja)%area)
2446 dnorm = norm2(cts(ia)%normal)
2447 cts(ia)%normal = cts(ia)%normal/dnorm
2448
2450 cts(ia)%r_sphere = (cts(ia)%r_sphere*cts(ia)%area + cts(ja)%r_sphere*cts(ja)%area) / &
2451 (cts(ia)%area + cts(ja)%area)
2452
2454 cts(ia)%area = cts(ia)%area + cts(ja)%area
2455
2457 do ii = ja+1, nts
2458 cts(ii-1) = cts(ii)
2459 end do
2460 nts = nts - 1 ! updating number of tesserae
2461 band_iter = .false.
2462 exit loop_ia
2463
2464 end do loop_ja
2465 end do loop_ia
2466 end do
2467
2469 vol = m_zero
2470 do its = 1, nts
2471 prod = dot_product(cts(its)%point, cts(its)%normal)
2472 vol = vol + cts(its)%area * prod / m_three
2473 stot = stot + cts(its)%area
2474 end do
2475
2476 if (mpi_world%is_root()) then
2477 write(unit_pcminfo, '(A2)') '# '
2478 write(unit_pcminfo, '(A29,I4)') '# Total number of tesserae = ' , nts
2479 write(unit_pcminfo, '(A30,F12.6)') '# Cavity surface area (A^2) = ', stot
2480 write(unit_pcminfo, '(A24,F12.6)') '# Cavity volume (A^3) = ' , vol
2481 write(unit_pcminfo, '(A2)') '# '
2482 end if
2483
2485 cts(:)%area = cts(:)%area*p_ang*p_ang
2486 cts(:)%point(1) = cts(:)%point(1)*p_ang
2487 cts(:)%point(2) = cts(:)%point(2)*p_ang
2488 cts(:)%point(3) = cts(:)%point(3)*p_ang
2489 cts(:)%r_sphere = cts(:)%r_sphere*p_ang
2490
2491 sfe(:)%x=sfe(:)%x*p_ang
2492 sfe(:)%y=sfe(:)%y*p_ang
2493 sfe(:)%z=sfe(:)%z*p_ang
2494 sfe(:)%r=sfe(:)%r*p_ang
2495
2496 pop_sub(cav_gen)
2497 end subroutine cav_gen
2498
2499 ! -----------------------------------------------------------------------------
2500
2503 subroutine subtessera(sfe, ns, nesf, nv, pts, ccc, pp, pp1, area)
2504 type(pcm_sphere_t), intent(in) :: sfe(:)
2505 integer, intent(in) :: ns
2506 integer, intent(in) :: nesf
2507 integer, intent(inout) :: nv
2508 real(real64), intent(inout) :: pts(:,:)
2509 real(real64), intent(out) :: ccc(:,:)
2510 real(real64), intent(out) :: pp(:)
2511 real(real64), intent(out) :: pp1(:)
2512 real(real64), intent(out) :: area
2513
2514 real(real64), parameter :: TOL = -1.0e-10_real64
2515 integer, parameter :: DIM_TEN = 10
2516
2517 integer :: intsph(1:DIM_TEN)
2518 integer :: nsfe1
2519 integer :: na
2520 integer :: icop
2521 integer :: ll
2522 integer :: iv1
2523 integer :: iv2
2524 integer :: ii
2525 integer :: icut
2526 integer :: jj
2527
2528 real(real64) :: p1(1:PCM_DIM_SPACE)
2529 real(real64) :: p2(1:PCM_DIM_SPACE)
2530 real(real64) :: p3(1:PCM_DIM_SPACE)
2531 real(real64) :: p4(1:PCM_DIM_SPACE)
2532 real(real64) :: point(1:PCM_DIM_SPACE)
2533 real(real64) :: pscr(1:PCM_DIM_SPACE,1:DIM_TEN)
2534 real(real64) :: cccp(1:PCM_DIM_SPACE,1:DIM_TEN)
2535 real(real64) :: pointl(1:PCM_DIM_SPACE,1:DIM_TEN)
2536 real(real64) :: diff(1:PCM_DIM_SPACE)
2537
2538 integer :: ind(1:DIM_TEN)
2539 integer :: ltyp(1:DIM_TEN)
2540 integer :: intscr(1:DIM_TEN)
2541
2542 real(real64) :: delr
2543 real(real64) :: delr2
2544 real(real64) :: rc
2545 real(real64) :: rc2
2546 real(real64) :: dnorm
2547 real(real64) :: dist
2548 real(real64) :: de2
2549
2550 push_sub(subtessera)
2551
2552 p1 = m_zero
2553 p2 = m_zero
2554 p3 = m_zero
2555 p4 = m_zero
2556 point = m_zero
2557 pscr = m_zero
2558 cccp = m_zero
2559 pointl = m_zero
2560 diff = m_zero
2561 area = m_zero
2562
2563 do jj = 1, 3
2564 ccc(1,jj) = sfe(ns)%x
2565 ccc(2,jj) = sfe(ns)%y
2566 ccc(3,jj) = sfe(ns)%z
2567 end do
2568
2569 intsph = ns
2570 do nsfe1 = 1, nesf
2571 if (nsfe1 == ns) cycle
2572 do jj = 1, nv
2573 intscr(jj) = intsph(jj)
2574 pscr(:,jj) = pts(:,jj)
2575 cccp(:,jj) = ccc(:,jj)
2576 end do
2577
2578 icop = 0
2579 ind = 0
2580 ltyp = 0
2581
2582 do ii = 1, nv
2583 delr2 = (pts(1,ii) - sfe(nsfe1)%x)**2 + (pts(2,ii) - sfe(nsfe1)%y)**2 + (pts(3,ii) - sfe(nsfe1)%z)**2
2584 delr = sqrt(delr2)
2585 if (delr < sfe(nsfe1)%r) then
2586 ind(ii) = 1
2587 icop = icop + 1
2588 end if
2589 end do
2590
2591 if (icop == nv) then
2592 pop_sub(subtessera)
2593 return
2594 end if
2595
2596 do ll = 1, nv
2597 iv1 = ll
2598 iv2 = ll+1
2599 if (ll == nv) iv2 = 1
2600 if ((ind(iv1) == 1) .and. (ind(iv2) == 1)) then
2601 ltyp(ll) = 0
2602 else if ((ind(iv1) == 0) .and. (ind(iv2) == 1)) then
2603 ltyp(ll) = 1
2604 else if ((ind(iv1) == 1) .and. (ind(iv2) == 0)) then
2605 ltyp(ll) = 2
2606 else if ((ind(iv1) == 0) .and. (ind(iv2) == 0)) then
2607 ltyp(ll) = 4
2608 diff = ccc(:,ll) - pts(:,ll)
2609 rc2 = dot_product(diff, diff)
2610 rc = sqrt(rc2)
2611
2612 do ii = 1, 11
2613 point = pts(:,iv1) + ii * (pts(:,iv2) - pts(:,iv1)) / 11
2614 point = point - ccc(:,ll)
2615 dnorm = norm2(point)
2616 point = point * rc / dnorm + ccc(:,ll)
2617
2618 dist = sqrt((point(1) - sfe(nsfe1)%x)**2 + (point(2) - sfe(nsfe1)%y)**2 + (point(3) - sfe(nsfe1)%z)**2)
2619
2620 if ((dist - sfe(nsfe1)%r) < tol) then
2621 ltyp(ll) = 3
2622 pointl(:, ll) = point
2623 exit
2624 end if
2625
2626 end do
2627 end if
2628 end do
2629
2630 icut = 0
2631 do ll = 1, nv
2632 if ((ltyp(ll) == 1) .or. (ltyp(ll) == 2)) icut = icut + 1
2633 if (ltyp(ll) == 3) icut = icut + 2
2634 end do
2635 icut = icut / 2
2636 if (icut > 1) then
2637 pop_sub(subtessera)
2638 return
2639 end if
2640
2641 na = 1
2642 do ll = 1, nv
2643
2644 if (ltyp(ll) == 0) cycle
2645 iv1 = ll
2646 iv2 = ll + 1
2647 if (ll == nv) iv2 = 1
2648
2649 if (ltyp(ll) == 1) then
2650 pts(:,na) = pscr(:,iv1)
2651 ccc(:,na) = cccp(:,iv1)
2652 intsph(na) = intscr(iv1)
2653 na = na + 1
2654 p1 = pscr(:,iv1)
2655 p2 = pscr(:,iv2)
2656 p3 = cccp(:,iv1)
2657
2658 call inter(sfe, p1, p2, p3, p4, nsfe1, 0)
2659 pts(:,na) = p4
2660
2661 de2 = (sfe(nsfe1)%x - sfe(ns)%x)**2 + ( sfe(nsfe1)%y - sfe(ns)%y)**2 + &
2662 (sfe(nsfe1)%z - sfe(ns)%z)**2
2663
2664 ccc(1,na) = sfe(ns)%x + ( sfe(nsfe1)%x - sfe(ns)%x)* &
2665 (sfe(ns)%r**2 - sfe(nsfe1)%r**2 + de2) / (m_two*de2)
2666
2667 ccc(2,na) = sfe(ns)%y + ( sfe(nsfe1)%y - sfe(ns)%y)* &
2668 (sfe(ns)%r**2 - sfe(nsfe1)%r**2 + de2) / (m_two*de2)
2669
2670 ccc(3,na) = sfe(ns)%z + ( sfe(nsfe1)%z - sfe(ns)%z)* &
2671 (sfe(ns)%r**2 - sfe(nsfe1)%r**2 + de2) / (m_two*de2)
2672
2673 intsph(na) = nsfe1
2674 na = na + 1
2675 end if
2676
2677 if (ltyp(ll) == 2) then
2678 p1 = pscr(:,iv1)
2679 p2 = pscr(:,iv2)
2680 p3 = cccp(:,iv1)
2681
2682 call inter(sfe, p1, p2, p3, p4, nsfe1, 1)
2683 pts(:,na) = p4
2684 ccc(:,na) = cccp(:,iv1)
2685 intsph(na) = intscr(iv1)
2686 na = na + 1
2687 end if
2688
2689 if (ltyp(ll) == 3) then
2690 pts(:,na) = pscr(:,iv1)
2691 ccc(:,na) = cccp(:,iv1)
2692 intsph(na) = intscr(iv1)
2693 na = na + 1
2694 p1 = pscr(:,iv1)
2695 p2 = pointl(:,ll)
2696 p3 = cccp(:,iv1)
2697
2698 call inter(sfe, p1, p2, p3, p4, nsfe1, 0)
2699 pts(:,na) = p4
2700
2701 de2 = (sfe(nsfe1)%x - sfe(ns)%x)**2 + (sfe(nsfe1)%y - sfe(ns)%y)**2 + (sfe(nsfe1)%z - sfe(ns)%z)**2
2702
2703 ccc(1,na) = sfe(ns)%x + (sfe(nsfe1)%x - sfe(ns)%x) * (sfe(ns)%r**2 - sfe(nsfe1)%r**2 + de2) / (m_two*de2)
2704
2705 ccc(2,na) = sfe(ns)%y + (sfe(nsfe1)%y - sfe(ns)%y) * (sfe(ns)%r**2 - sfe(nsfe1)%r**2 + de2) / (m_two*de2)
2706
2707 ccc(3,na) = sfe(ns)%z + (sfe(nsfe1)%z - sfe(ns)%z) * (sfe(ns)%r**2 - sfe(nsfe1)%r**2 + de2) / (m_two*de2)
2708
2709 intsph(na) = nsfe1
2710 na = na + 1
2711 p1 = pointl(:,ll)
2712 p2 = pscr(:,iv2)
2713 p3 = cccp(:,iv1)
2714
2715 call inter(sfe, p1, p2, p3, p4, nsfe1, 1)
2716 pts(:,na) = p4
2717 ccc(:,na) = cccp(:,iv1)
2718 intsph(na) = intscr(iv1)
2719 na = na + 1
2720 end if
2721
2722 if (ltyp(ll) == 4) then
2723 pts(:,na) = pscr(:,iv1)
2724 ccc(:,na) = cccp(:,iv1)
2725 intsph(na) = intscr(iv1)
2726 na = na + 1
2727 end if
2728 end do
2729
2730 nv = na - 1
2731 if (nv > 10) then
2732 message(1) = "Too many vertices on the tessera"
2733 call messages_fatal(1)
2734 end if
2735 end do
2736
2737 call gaubon(sfe, nv, ns, pts, ccc, pp, pp1, area, intsph)
2738
2739 pop_sub(subtessera)
2740 end subroutine subtessera
2741
2742 ! -----------------------------------------------------------------------------
2743
2747 subroutine inter(sfe, p1, p2, p3, p4, ns, ia)
2748 type(pcm_sphere_t), intent(in) :: sfe(:)
2749 real(real64), intent(in) :: p1(1:PCM_DIM_SPACE)
2750 real(real64), intent(in) :: p2(1:PCM_DIM_SPACE)
2751 real(real64), intent(in) :: p3(1:PCM_DIM_SPACE)
2752 real(real64), intent(out) :: p4(1:PCM_DIM_SPACE)
2753 integer, intent(in) :: ns
2754 integer, intent(in) :: ia
2755
2756 real(real64), parameter :: TOL = 1.0e-08_real64
2757
2758 integer :: m_iter
2759 real(real64) :: r
2760 real(real64) :: alpha
2761 real(real64) :: delta
2762 real(real64) :: dnorm
2763 real(real64) :: diff
2764 real(real64) :: diff_vec(1:PCM_DIM_SPACE)
2765 logical :: band_iter
2766
2767 push_sub(inter)
2768
2769 diff_vec = p1 - p3
2770 r = norm2(diff_vec)
2771
2772 alpha = m_half
2773 delta = m_zero
2774 m_iter = 1
2775
2776 band_iter = .false.
2777 do while(.not.(band_iter))
2778 if (m_iter > 1000) then
2779 message(1) = "Too many iterations inside subrotuine inter"
2780 call messages_fatal(1)
2781 end if
2782
2783 band_iter = .true.
2784
2785 alpha = alpha + delta
2786
2787 p4 = p1 + alpha*(p2 - p1) - p3
2788 dnorm = norm2(p4)
2789 p4 = p4*r/dnorm + p3
2790 diff = (p4(1) - sfe(ns)%x)**2 + (p4(2) - sfe(ns)%y)**2 + (p4(3) - sfe(ns)%z)**2
2791 diff = sqrt(diff) - sfe(ns)%r
2792
2793 if (abs(diff) < tol) then
2794 pop_sub(inter)
2795 return
2796 end if
2797
2798 if (ia == 0) then
2799 if (diff > m_zero) delta = m_one/(m_two**(m_iter + 1))
2800 if (diff < m_zero) delta = -m_one/(m_two**(m_iter + 1))
2801 m_iter = m_iter + 1
2802 band_iter = .false.
2803 end if
2804
2805 if (ia == 1) then
2806 if (diff > m_zero) delta = -m_one/(m_two**(m_iter + 1))
2807 if (diff < m_zero) delta = m_one/(m_two**(m_iter + 1))
2808 m_iter = m_iter + 1
2809 band_iter = .false.
2810 end if
2811 end do
2812
2813 pop_sub(inter)
2814 end subroutine inter
2815
2816 ! -----------------------------------------------------------------------------
2817
2824 subroutine gaubon(sfe, nv, ns, pts, ccc, pp, pp1, area, intsph)
2825 type(pcm_sphere_t), intent(in) :: sfe(:)
2826 real(real64), intent(in) :: pts(:,:)
2827 real(real64), intent(in) :: ccc(:,:)
2828 real(real64), intent(inout) :: pp(:)
2829 real(real64), intent(inout) :: pp1(:)
2830 integer, intent(in) :: intsph(:)
2831 real(real64), intent(out) :: area
2832 integer, intent(in) :: nv
2833 integer, intent(in) :: ns
2834
2835 real(real64) :: p1(1:PCM_DIM_SPACE), p2(1:PCM_DIM_SPACE), p3(1:PCM_DIM_SPACE)
2836 real(real64) :: u1(1:PCM_DIM_SPACE), u2(1:PCM_DIM_SPACE)
2837 real(real64) :: point_1(1:PCM_DIM_SPACE), point_2(1:PCM_DIM_SPACE)
2838 real(real64) :: tpi, sum1, dnorm, dnorm1, dnorm2
2839 real(real64) :: cosphin, phin, costn, sum2, betan
2840 integer :: nsfe1, ia, nn, n0, n1, n2
2841
2842 push_sub(gaubon)
2843
2844 point_1 = m_zero
2845 point_2 = m_zero
2846 p1 = m_zero
2847 p2 = m_zero
2848 p3 = m_zero
2849 u1 = m_zero
2850 u2 = m_zero
2851
2852 tpi = m_two*m_pi
2853 sum1 = m_zero
2854 do nn = 1, nv
2855 point_1 = pts(:,nn) - ccc(:,nn)
2856 if (nn < nv) then
2857 point_2 = pts(:,nn+1) - ccc(:,nn)
2858 else
2859 point_2 = pts(:,1) - ccc(:,nn)
2860 end if
2861
2862 dnorm1 = norm2(point_1)
2863 dnorm2 = norm2(point_2)
2864 cosphin = dot_product(point_1, point_2) / (dnorm1*dnorm2)
2865
2866 if (cosphin > m_one) cosphin = m_one
2867 if (cosphin < -m_one) cosphin = -m_one
2868
2869 phin = acos(cosphin)
2870 nsfe1 = intsph(nn)
2871
2872 point_1(1) = sfe(nsfe1)%x - sfe(ns)%x
2873 point_1(2) = sfe(nsfe1)%y - sfe(ns)%y
2874 point_1(3) = sfe(nsfe1)%z - sfe(ns)%z
2875
2876 dnorm1 = norm2(point_1)
2877
2878 if (abs(dnorm1) <= m_epsilon) dnorm1 = m_one
2879
2880 point_2(1) = pts(1,nn) - sfe(ns)%x
2881 point_2(2) = pts(2,nn) - sfe(ns)%y
2882 point_2(3) = pts(3,nn) - sfe(ns)%z
2883
2884 dnorm2 = norm2(point_2)
2885
2886 costn = dot_product(point_1, point_2) / (dnorm1 * dnorm2)
2887 sum1 = sum1 + phin * costn
2888 end do
2889
2890 sum2 = m_zero
2892 do nn = 1, nv
2893 p1 = m_zero
2894 p2 = m_zero
2895 p3 = m_zero
2896
2897 n1 = nn
2898 if (nn > 1) n0 = nn - 1
2899 if (nn == 1) n0 = nv
2900 if (nn < nv) n2 = nn + 1
2901 if (nn == nv) n2 = 1
2902
2903 p1 = pts(:,n1) - ccc(:,n0)
2904 p2 = pts(:,n0) - ccc(:,n0)
2905 call vecp(p1, p2, p3, dnorm)
2906 p2 = p3
2907
2908 call vecp(p1, p2, p3, dnorm)
2909 u1 = p3/dnorm
2910
2911 p1 = pts(:,n1) - ccc(:,n1)
2912 p2 = pts(:,n2) - ccc(:,n1)
2913 call vecp(p1, p2, p3, dnorm)
2914 p2 = p3
2915
2916 call vecp(p1, p2, p3, dnorm)
2917 u2 = p3/dnorm
2918
2919 betan = acos(dot_product(u1, u2))
2920 sum2 = sum2 + (m_pi - betan)
2921 end do
2922
2924 area = sfe(ns)%r*sfe(ns)%r*(tpi + sum1 - sum2)
2925
2927 pp = m_zero
2928
2929 do ia = 1, nv
2930 pp(1) = pp(1) + (pts(1,ia) - sfe(ns)%x)
2931 pp(2) = pp(2) + (pts(2,ia) - sfe(ns)%y)
2932 pp(3) = pp(3) + (pts(3,ia) - sfe(ns)%z)
2933 end do
2934
2935 dnorm = norm2(pp)
2936
2937 pp(1) = sfe(ns)%x + pp(1) * sfe(ns)%r / dnorm
2938 pp(2) = sfe(ns)%y + pp(2) * sfe(ns)%r / dnorm
2939 pp(3) = sfe(ns)%z + pp(3) * sfe(ns)%r / dnorm
2940
2942 pp1(1) = (pp(1) - sfe(ns)%x) / sfe(ns)%r
2943 pp1(2) = (pp(2) - sfe(ns)%y) / sfe(ns)%r
2944 pp1(3) = (pp(3) - sfe(ns)%z) / sfe(ns)%r
2945
2947 if (area < m_zero) area = m_zero
2948
2949 pop_sub(gaubon)
2950 end subroutine gaubon
2951
2952 ! -----------------------------------------------------------------------------
2953
2955 subroutine vecp(p1, p2, p3, dnorm)
2956 real(real64), intent(in) :: P1(:)
2957 real(real64), intent(in) :: P2(:)
2958 real(real64), intent(out) :: P3(:)
2959 real(real64), intent(out) :: dnorm
2960
2961 p3 = m_zero
2962 p3(1) = p1(2)*p2(3) - p1(3)*p2(2)
2963 p3(2) = p1(3)*p2(1) - p1(1)*p2(3)
2964 p3(3) = p1(1)*p2(2) - p1(2)*p2(1)
2965
2966 dnorm = norm2(p3)
2967
2968 end subroutine vecp
2969
2970 subroutine pcm_end(pcm)
2971 type(pcm_t), intent(inout) :: pcm
2972
2973 integer :: asc_unit_test, asc_unit_test_sol, asc_unit_test_e, asc_unit_test_n, asc_unit_test_ext
2974 integer :: ia
2975
2976 push_sub(pcm_end)
2977
2978 if (pcm%solute .and. pcm%localf) then
2979 asc_unit_test = io_open(pcm_dir//'ASC.dat', pcm%namespace, action='write')
2980 asc_unit_test_sol = io_open(pcm_dir//'ASC_sol.dat', pcm%namespace, action='write')
2981 asc_unit_test_e = io_open(pcm_dir//'ASC_e.dat', pcm%namespace, action='write')
2982 asc_unit_test_n = io_open(pcm_dir//'ASC_n.dat', pcm%namespace, action='write')
2983 asc_unit_test_ext = io_open(pcm_dir//'ASC_ext.dat', pcm%namespace, action='write')
2984 do ia = 1, pcm%n_tesserae
2985 write(asc_unit_test,*) pcm%tess(ia)%point, pcm%q_e(ia) + pcm%q_n(ia) + pcm%q_ext(ia), ia
2986 write(asc_unit_test_sol,*) pcm%tess(ia)%point, pcm%q_e(ia) + pcm%q_n(ia), ia
2987 write(asc_unit_test_e,*) pcm%tess(ia)%point, pcm%q_e(ia), ia
2988 write(asc_unit_test_n,*) pcm%tess(ia)%point, pcm%q_n(ia), ia
2989 write(asc_unit_test_ext,*) pcm%tess(ia)%point, pcm%q_ext(ia), ia
2990 end do
2991 call io_close(asc_unit_test)
2992 call io_close(asc_unit_test_sol)
2993 call io_close(asc_unit_test_e)
2994 call io_close(asc_unit_test_n)
2995 call io_close(asc_unit_test_ext)
2996
2997 else if (pcm%solute .and. .not. pcm%localf) then
2998 asc_unit_test_sol = io_open(pcm_dir//'ASC_sol.dat', pcm%namespace, action='write')
2999 asc_unit_test_e = io_open(pcm_dir//'ASC_e.dat', pcm%namespace, action='write')
3000 asc_unit_test_n = io_open(pcm_dir//'ASC_n.dat', pcm%namespace, action='write')
3001 do ia = 1, pcm%n_tesserae
3002 write(asc_unit_test_sol,*) pcm%tess(ia)%point, pcm%q_e(ia) + pcm%q_n(ia), ia
3003 write(asc_unit_test_e,*) pcm%tess(ia)%point, pcm%q_e(ia), ia
3004 write(asc_unit_test_n,*) pcm%tess(ia)%point, pcm%q_n(ia), ia
3005 end do
3006 call io_close(asc_unit_test_sol)
3007 call io_close(asc_unit_test_e)
3008 call io_close(asc_unit_test_n)
3009
3010 else if (.not. pcm%solute .and. pcm%localf) then
3011 asc_unit_test_ext = io_open(pcm_dir//'ASC_ext.dat', pcm%namespace, action='write')
3012 do ia = 1, pcm%n_tesserae
3013 write(asc_unit_test_ext,*) pcm%tess(ia)%point, pcm%q_ext(ia), ia
3014 end do
3015 call io_close(asc_unit_test_ext)
3016 end if
3017
3018 safe_deallocate_a(pcm%spheres)
3019 safe_deallocate_a(pcm%tess)
3020 safe_deallocate_a(pcm%matrix)
3021 if (.not. is_close(pcm%epsilon_infty, pcm%epsilon_0)) then
3022 safe_deallocate_a(pcm%matrix_d)
3023 end if
3024 safe_deallocate_a(pcm%q_e)
3025 safe_deallocate_a(pcm%q_e_in)
3026 safe_deallocate_a(pcm%q_n)
3027 safe_deallocate_a(pcm%v_e)
3028 safe_deallocate_a(pcm%v_n)
3029 safe_deallocate_a(pcm%v_e_rs)
3030 safe_deallocate_a(pcm%v_n_rs)
3031 if (pcm%localf) then
3032 safe_deallocate_a(pcm%matrix_lf)
3033 if (.not. is_close(pcm%epsilon_infty, pcm%epsilon_0)) then
3034 safe_deallocate_a(pcm%matrix_lf_d)
3035 end if
3036
3037 safe_deallocate_a(pcm%q_ext)
3038 safe_deallocate_a(pcm%q_ext_in)
3039 safe_deallocate_a(pcm%v_ext)
3040 safe_deallocate_a(pcm%v_ext_rs)
3041 if (pcm%kick_is_present) then
3042 safe_deallocate_a(pcm%q_kick)
3043 safe_deallocate_a(pcm%v_kick)
3044 safe_deallocate_a(pcm%v_kick_rs)
3045 end if
3046 end if
3047
3048 if (pcm%calc_method == pcm_calc_poisson) then
3049 safe_deallocate_a( pcm%rho_n)
3050 safe_deallocate_a( pcm%rho_e)
3051 if (pcm%localf) then
3052 safe_deallocate_a( pcm%rho_ext)
3053 if (pcm%kick_is_present) then
3054 safe_deallocate_a( pcm%rho_kick)
3055 end if
3056 end if
3057 end if
3058
3059 if (pcm%tdlevel == pcm_td_eom) call pcm_eom_end()
3060
3061 if (mpi_world%is_root()) call io_close(pcm%info_unit)
3062
3063 pop_sub(pcm_end)
3064 end subroutine pcm_end
3065
3066 ! -----------------------------------------------------------------------------
3068 logical function pcm_update(this) result(update)
3069 type(pcm_t), intent(inout) :: this
3070
3071 this%iter = this%iter + 1
3072 update = this%iter <= 6 .or. mod(this%iter, this%update_iter) == 0
3073
3074 end function pcm_update
3075
3076 ! -----------------------------------------------------------------------------
3078 real(real64) function pcm_get_vdw_radius(species, pcm_vdw_type, namespace) result(vdw_r)
3079 class(species_t), intent(in) :: species
3080 integer, intent(in) :: pcm_vdw_type
3081 type(namespace_t), intent(in) :: namespace
3082
3083 integer :: ia
3084 integer, parameter :: upto_xe = 54
3085 real(real64), save :: vdw_radii(1:upto_xe)
3086 ! by Stefan Grimme in J. Comput. Chem. 27: 1787-1799, 2006
3087 ! except for C, N and O, reported in J. Chem. Phys. 120, 3893 (2004).
3088 data (vdw_radii(ia), ia=1, upto_xe) / &
3089 !H He
3090 1.001_real64, 1.012_real64, &
3091 !Li Be B C N O F Ne
3092 0.825_real64, 1.408_real64, 1.485_real64, 2.000_real64, 1.583_real64, 1.500_real64, 1.287_real64, 1.243_real64, &
3093 !Na Mg Al Si P S Cl Ar
3094 1.144_real64, 1.364_real64, 1.639_real64, 1.716_real64, 1.705_real64, 1.683_real64, 1.639_real64, 1.595_real64, &
3095 !K Ca
3096 1.485_real64, 1.474_real64, &
3098 1.562_real64, 1.562_real64, &
3099 1.562_real64, 1.562_real64, &
3100 1.562_real64, 1.562_real64, &
3101 1.562_real64, 1.562_real64, &
3102 1.562_real64, 1.562_real64, &
3103 !Ga Ge As Se Br Kr
3104 1.650_real64, 1.727_real64, 1.760_real64, 1.771_real64, 1.749_real64, 1.727_real64, &
3105 !Rb Sr !> Y -- Cd <!
3106 1.628_real64, 1.606_real64, 1.639_real64, 1.639_real64, &
3107 1.639_real64, 1.639_real64, &
3108 1.639_real64, 1.639_real64, &
3109 1.639_real64, 1.639_real64, &
3110 1.639_real64, 1.639_real64, &
3111 !In Sn Sb Te I Xe
3112 2.672_real64, 1.804_real64, 1.881_real64, 1.892_real64, 1.892_real64, 1.881_real64 /
3113
3114 select case (pcm_vdw_type)
3115 case (pcm_vdw_optimized)
3116 if (species%get_z() > upto_xe) then
3117 write(message(1),'(a,a)') "The van der Waals radius is missing for element ", trim(species%get_label())
3118 write(message(2),'(a)') "Use PCMVdWRadii = pcm_vdw_species, for other vdw radii values"
3119 call messages_fatal(2, namespace=namespace)
3120 end if
3121 ia = int(species%get_z())
3122 vdw_r = vdw_radii(ia)*p_ang
3123
3124 case (pcm_vdw_species)
3125 vdw_r = species%get_vdw_radius()
3126 if (vdw_r < m_zero) then
3127 call messages_write('The default vdW radius for species '//trim(species%get_label())//':')
3128 call messages_write(' is not defined. ')
3129 call messages_write(' Add a positive vdW radius value in %Species block. ')
3130 call messages_fatal(namespace=namespace)
3131 end if
3132 end select
3133
3134 end function pcm_get_vdw_radius
3135
3136
3137 ! -----------------------------------------------------------------------------
3139 subroutine pcm_dipole(mu_pcm, q_pcm, tess, n_tess)
3140 real(real64), intent(out) :: mu_pcm(:)
3141 real(real64), intent(in) :: q_pcm(:)
3142 integer, intent(in) :: n_tess
3143 type(pcm_tessera_t), intent(in) :: tess(:)
3144
3145 integer :: ia
3146
3147 push_sub(pcm_dipole)
3148
3149 mu_pcm = m_zero
3150 do ia = 1, n_tess
3151 mu_pcm = mu_pcm + q_pcm(ia) * tess(ia)%point
3152 end do
3153
3154 pop_sub(pcm_dipole)
3155 end subroutine pcm_dipole
3156
3157 ! -----------------------------------------------------------------------------
3158 ! Driver function to evaluate eps(omega)
3159 subroutine pcm_eps(pcm, eps, omega)
3160 type(pcm_min_t), intent(in) :: pcm
3161 complex(real64), intent(out) :: eps
3162 real(real64), intent(in) :: omega
3163
3164 push_sub(pcm_eps)
3165
3166 if (pcm%tdlevel == pcm_td_eom) then
3167 if (pcm%which_eps == pcm_debye_model) then
3168 call pcm_eps_deb(eps, pcm%deb, omega)
3169 else if (pcm%which_eps == pcm_drude_model) then
3170 call pcm_eps_drl(eps, pcm%drl, omega)
3171 end if
3172 else if (pcm%tdlevel == pcm_td_neq) then
3173 eps = pcm%deb%eps_d
3174 else if (pcm%tdlevel == pcm_td_eq) then
3175 eps = pcm%deb%eps_0
3176 end if
3177
3178 pop_sub(pcm_eps)
3179 end subroutine pcm_eps
3180
3181 ! -----------------------------------------------------------------------------
3182
3183 subroutine pcm_min_input_parsing_for_spectrum(pcm, namespace)
3184 type(pcm_min_t), intent(out) :: pcm
3185 type(namespace_t), intent(in) :: namespace
3186
3187 push_sub(pcm_min_input_parsing_for_spectrum)
3188
3189 ! re-parsing PCM keywords
3190 call parse_variable(namespace, 'PCMCalculation', .false., pcm%run_pcm)
3191 call messages_print_with_emphasis(msg='PCM', namespace=namespace)
3192 call parse_variable(namespace, 'PCMLocalField', .false., pcm%localf)
3193 call messages_print_var_value("PCMLocalField", pcm%localf, namespace=namespace)
3194 if (pcm%localf) then
3195 call messages_experimental("PCM local field effects in the optical spectrum", namespace=namespace)
3196 call messages_write('Beware of possible numerical errors in the optical spectrum due to PCM local field effects,')
3197 call messages_new_line()
3198 call messages_write('particularly, when static and high-frequency values of the dielectric functions are large')
3199 call messages_write(' (>~10 in units of the vacuum permittivity \epsilon_0).')
3200 call messages_new_line()
3201 call messages_write('However, PCM local field effects in the optical spectrum work well for polar or non-polar solvents')
3202 call messages_new_line()
3203 call messages_write('in the nonequilibrium or equation-of-motion TD-PCM propagation schemes.')
3204 call messages_warning(namespace=namespace)
3205 end if
3206 call parse_variable(namespace, 'PCMTDLevel' , pcm_td_eq, pcm%tdlevel)
3207 call messages_print_var_value("PCMTDLevel", pcm%tdlevel, namespace=namespace)
3208
3209 ! reading dielectric function model parameters
3210 call parse_variable(namespace, 'PCMStaticEpsilon' , m_one, pcm%deb%eps_0)
3211 call messages_print_var_value("PCMStaticEpsilon", pcm%deb%eps_0, namespace=namespace)
3212 if (pcm%tdlevel == pcm_td_eom) then
3213 call parse_variable(namespace, 'PCMEpsilonModel', pcm_debye_model, pcm%which_eps)
3214 call messages_print_var_value("PCMEpsilonModel", pcm%which_eps, namespace=namespace)
3215 if (pcm%which_eps == pcm_debye_model) then
3216 call parse_variable(namespace, 'PCMDynamicEpsilon', pcm%deb%eps_0, pcm%deb%eps_d)
3217 call messages_print_var_value("PCMDynamicEpsilon", pcm%deb%eps_d, namespace=namespace)
3218 call parse_variable(namespace, 'PCMDebyeRelaxTime', m_zero, pcm%deb%tau)
3219 call messages_print_var_value("PCMDebyeRelaxTime", pcm%deb%tau, namespace=namespace)
3220 else if (pcm%which_eps == pcm_drude_model) then
3221 call parse_variable(namespace, 'PCMDrudeLOmega', sqrt(m_one/(pcm%deb%eps_0-m_one)), pcm%drl%w0)
3222 call messages_print_var_value("PCMDrudeLOmega", pcm%drl%w0, namespace=namespace)
3223 call parse_variable(namespace, 'PCMDrudeLDamping', m_zero, pcm%drl%gm)
3224 call messages_print_var_value("PCMDrudeLDamping", pcm%drl%gm, namespace=namespace)
3225 end if
3226 else if (pcm%tdlevel == pcm_td_neq) then
3227 call parse_variable(namespace, 'PCMDynamicEpsilon', pcm%deb%eps_0, pcm%deb%eps_d)
3228 call messages_print_var_value("PCMDynamicEpsilon", pcm%deb%eps_d, namespace=namespace)
3229 end if
3230
3231 pop_sub(pcm_min_input_parsing_for_spectrum)
3232 end subroutine pcm_min_input_parsing_for_spectrum
3233
3234 ! -----------------------------------------------------------------------------
3235 ! Debye dielectric function
3236 subroutine pcm_eps_deb(eps, deb, omega)
3237 complex(real64), intent(out) :: eps
3238 type(debye_param_t), intent(in) :: deb
3239 real(real64), intent(in) :: omega
3240
3241 push_sub(pcm_eps_deb)
3242
3243 eps = deb%eps_d + (deb%eps_0 - deb%eps_d)/(m_one + (omega*deb%tau)**2) + &
3244 m_zi*omega*deb%tau*(deb%eps_0 - deb%eps_d)/(m_one + (omega*deb%tau)**2)
3245
3246 pop_sub(pcm_eps_deb)
3247 end subroutine pcm_eps_deb
3248
3249 ! -----------------------------------------------------------------------------
3250 ! Drude-Lorentz dielectric function
3251 subroutine pcm_eps_drl(eps, drl, omega)
3252 complex(real64), intent(out) :: eps
3253 type(drude_param_t), intent(in) :: drl
3254 real(real64), intent(in) :: omega
3255
3256 push_sub(pcm_eps_drl)
3257
3258 eps = m_one + (drl%w0**2 - omega**2)*drl%aa/((drl%w0**2 - omega**2)**2 + (omega*drl%gm)**2) + &
3259 m_zi*omega*drl%gm*drl%aa/((drl%w0**2 - omega**2)**2 + (omega*drl%gm)**2)
3260
3261 pop_sub(pcm_eps_drl)
3262 end subroutine pcm_eps_drl
3263
3264end module pcm_oct_m
3265
3266!! Local Variables:
3267!! mode: f90
3268!! coding: utf-8
3269!! End:
info
subroutine info()
Definition: em_resp.F90:1093
lalg_basic_oct_m::lalg_gemv
Definition: lalg_basic.F90:217
math_oct_m::is_close
Definition: math.F90:171
mesh_function_oct_m::dmf_integrate
Definition: mesh_function.F90:207
mesh_interpolation_oct_m::mesh_interpolation_evaluate
Definition: mesh_interpolation.F90:147
messages_oct_m::messages_print_var_option
Definition: messages.F90:199
messages_oct_m::messages_print_var_value
Prints out to iunit a message in the form: ["InputVariable" = value] where "InputVariable" is given b...
Definition: messages.F90:182
messages_oct_m::messages_write
Definition: messages.F90:190
parser_oct_m::parse_variable
Definition: parser.F90:253
unit_oct_m::sqrt
Definition: unit.F90:176
unit_oct_m::units_from_atomic
Definition: unit.F90:172
varinfo_oct_m::varinfo_valid_option
Definition: varinfo.F90:134
acos
double acos(double __x) __attribute__((__nothrow__
exp
double exp(double __x) __attribute__((__nothrow__
sin
double sin(double __x) __attribute__((__nothrow__
sqrt
double sqrt(double __x) __attribute__((__nothrow__
cos
double cos(double __x) __attribute__((__nothrow__
floor
double floor(double __x) __attribute__((__nothrow__
box_minimum_oct_m
Definition: box_minimum.F90:117
debug_oct_m
Definition: debug.F90:116
debug_oct_m::debug
type(debug_t), save, public debug
Definition: debug.F90:158
global_oct_m
Definition: global.F90:116
global_oct_m::m_two
real(real64), parameter, public m_two
Definition: global.F90:202
global_oct_m::m_zero
real(real64), parameter, public m_zero
Definition: global.F90:200
global_oct_m::m_four
real(real64), parameter, public m_four
Definition: global.F90:204
global_oct_m::p_a_b
real(real64), parameter, public p_a_b
some physical constants
Definition: global.F90:237
global_oct_m::m_pi
real(real64), parameter, public m_pi
some mathematical constants
Definition: global.F90:198
global_oct_m::pcm_dir
character(len= *), parameter, public pcm_dir
Definition: global.F90:287
global_oct_m::m_epsilon
real(real64), parameter, public m_epsilon
Definition: global.F90:216
global_oct_m::m_half
real(real64), parameter, public m_half
Definition: global.F90:206
global_oct_m::m_one
real(real64), parameter, public m_one
Definition: global.F90:201
grid_oct_m
This module implements the underlying real-space grid.
Definition: grid.F90:119
index_oct_m
This module implements the index, used for the mesh points.
Definition: index.F90:124
io_function_oct_m
Definition: io_function.F90:116
io_function_oct_m::dio_function_output
subroutine, public dio_function_output(how, dir, fname, namespace, space, mesh, ff, unit, ierr, pos, atoms, grp, root)
Top-level IO routine for functions defined on the mesh.
Definition: io_function.F90:1684
io_function_oct_m::io_function_fill_how
integer(int64) function, public io_function_fill_how(where)
Use this function to quickly plot functions for debugging purposes: call dio_function_output(io_funct...
Definition: io_function.F90:825
io_oct_m
Definition: io.F90:116
io_oct_m::io_close
subroutine, public io_close(iunit, grp)
Definition: io.F90:467
io_oct_m::io_mkdir
subroutine, public io_mkdir(fname, namespace, parents)
Definition: io.F90:361
io_oct_m::io_open
integer function, public io_open(file, namespace, action, status, form, position, die, recl, grp)
Definition: io.F90:402
ions_oct_m
Definition: ions.F90:117
lalg_basic_oct_m
Definition: lalg_basic.F90:116
math_oct_m
This module is intended to contain "only mathematical" functions and procedures.
Definition: math.F90:117
mesh_function_oct_m
This module defines various routines, operating on mesh functions.
Definition: mesh_function.F90:118
mesh_interpolation_oct_m
Definition: mesh_interpolation.F90:116
mesh_interpolation_oct_m::mesh_interpolation_init
subroutine, public mesh_interpolation_init(this, mesh)
Definition: mesh_interpolation.F90:159
mesh_interpolation_oct_m::mesh_interpolation_end
subroutine, public mesh_interpolation_end(this)
Definition: mesh_interpolation.F90:173
mesh_oct_m
This module defines the meshes, which are used in Octopus.
Definition: mesh.F90:120
mesh_oct_m::mesh_local_index_from_coords
integer function, public mesh_local_index_from_coords(mesh, ix)
This function returns the local index of the point for a given vector of integer coordinates.
Definition: mesh.F90:938
mesh_oct_m::mesh_global_index_from_coords
integer(int64) function, public mesh_global_index_from_coords(mesh, ix)
This function returns the true global index of the point for a given vector of integer coordinates.
Definition: mesh.F90:919
mesh_oct_m::mesh_r
pure subroutine, public mesh_r(mesh, ip, rr, origin, coords)
return the distance to the origin for a given grid point
Definition: mesh.F90:342
messages_oct_m
Definition: messages.F90:117
messages_oct_m::messages_print_with_emphasis
subroutine, public messages_print_with_emphasis(msg, iunit, namespace)
Definition: messages.F90:898
messages_oct_m::msg
character(len=512), private msg
Definition: messages.F90:167
messages_oct_m::messages_warning
subroutine, public messages_warning(no_lines, all_nodes, namespace)
Definition: messages.F90:525
messages_oct_m::messages_new_line
subroutine, public messages_new_line()
Definition: messages.F90:1089
messages_oct_m::message
character(len=256), dimension(max_lines), public message
to be output by fatal, warning
Definition: messages.F90:162
messages_oct_m::messages_fatal
subroutine, public messages_fatal(no_lines, only_root_writes, namespace)
Definition: messages.F90:410
messages_oct_m::messages_input_error
subroutine, public messages_input_error(namespace, var, details, row, column)
Definition: messages.F90:691
messages_oct_m::messages_experimental
subroutine, public messages_experimental(name, namespace)
Definition: messages.F90:1040
messages_oct_m::messages_info
subroutine, public messages_info(no_lines, iunit, debug_only, stress, all_nodes, namespace)
Definition: messages.F90:594
mpi_oct_m
Definition: mpi.F90:116
mpi_oct_m::mpi_world
type(mpi_grp_t), public mpi_world
Definition: mpi.F90:272
namespace_oct_m
Definition: namespace.F90:105
par_vec_oct_m
Some general things and nomenclature:
Definition: par_vec.F90:173
parser_oct_m
Definition: parser.F90:116
pcm_eom_oct_m
Definition: pcm_eom.F90:115
pcm_eom_oct_m::pcm_external_plus_kick
integer, parameter, public pcm_external_plus_kick
Definition: pcm_eom.F90:174
pcm_eom_oct_m::pcm_electrons
integer, parameter, public pcm_electrons
Definition: pcm_eom.F90:174
pcm_eom_oct_m::pcm_kick
integer, parameter, public pcm_kick
Definition: pcm_eom.F90:174
pcm_eom_oct_m::pcm_nuclei
integer, parameter, public pcm_nuclei
Definition: pcm_eom.F90:174
pcm_eom_oct_m::pcm_charges_propagation
subroutine, public pcm_charges_propagation(q_t, pot_t, this_dt, this_cts_act, input_asc, this_eom, this_eps, namespace, this_deb, this_drl)
Driving subroutine for the Equation of Motion (EOM) propagation of the polarization charges within th...
Definition: pcm_eom.F90:235
pcm_eom_oct_m::pcm_external_potential
integer, parameter, public pcm_external_potential
Definition: pcm_eom.F90:174
pcm_eom_oct_m::pcm_drude_model
integer, parameter, public pcm_drude_model
Definition: pcm_eom.F90:170
pcm_eom_oct_m::pcm_eom_enough_initial
subroutine, public pcm_eom_enough_initial(not_yet_called)
Definition: pcm_eom.F90:940
pcm_eom_oct_m::pcm_debye_model
integer, parameter, public pcm_debye_model
Definition: pcm_eom.F90:170
pcm_oct_m
Definition: pcm.F90:116
pcm_oct_m::pcm_calc_pot_rs
subroutine, public pcm_calc_pot_rs(pcm, mesh, psolver, ions, v_h, v_ext, kick, time_present, kick_time)
Definition: pcm.F90:1217
pcm_oct_m::pcm_eps_deb
subroutine pcm_eps_deb(eps, deb, omega)
Definition: pcm.F90:3332
pcm_oct_m::pcm_update
logical function, public pcm_update(this)
Update pcm potential.
Definition: pcm.F90:3164
pcm_oct_m::pcm_pot_rs
subroutine, public pcm_pot_rs(pcm, v_pcm, q_pcm, rho, mesh, psolver)
Generates the potential 'v_pcm' in real-space.
Definition: pcm.F90:1885
pcm_oct_m::n_tess_sphere
integer, parameter n_tess_sphere
minimum number of tesserae per sphere
Definition: pcm.F90:288
pcm_oct_m::pcm_eps_drl
subroutine pcm_eps_drl(eps, drl, omega)
Definition: pcm.F90:3347
pcm_oct_m::pcm_nn_in_mesh
logical function pcm_nn_in_mesh(pcm, mesh)
Check wether the nearest neighbor requested are in the mesh or not.
Definition: pcm.F90:1688
pcm_oct_m::pcm_calc_direct
integer, parameter, public pcm_calc_direct
Definition: pcm.F90:280
pcm_oct_m::pcm_charge_density
subroutine, public pcm_charge_density(pcm, q_pcm, q_pcm_tot, mesh, rho)
Generates the polarization charge density smearing the charge with a gaussian distribution on the mes...
Definition: pcm.F90:1748
pcm_oct_m::pcm_get_vdw_radius
real(real64) function, public pcm_get_vdw_radius(species, pcm_vdw_type, namespace)
get the vdw radius
Definition: pcm.F90:3174
pcm_oct_m::pcm_matrix
subroutine pcm_matrix(eps, tess, n_tess, pcm_mat, localf)
Generates the PCM response matrix. J. Tomassi et al. Chem. Rev. 105, 2999 (2005).
Definition: pcm.F90:1994
pcm_oct_m::pcm_dipole
subroutine, public pcm_dipole(mu_pcm, q_pcm, tess, n_tess)
Computes the dipole moment mu_pcm due to a distribution of charges q_pcm.
Definition: pcm.F90:3235
pcm_oct_m::pcm_elect_energy
subroutine, public pcm_elect_energy(ions, pcm, E_int_ee, E_int_en, E_int_ne, E_int_nn, E_int_e_ext, E_int_n_ext)
Calculates the solute-solvent electrostatic interaction energy .
Definition: pcm.F90:1545
pcm_oct_m::s_mat_elem_i
real(real64) function s_mat_elem_i(tessi, tessj)
electrostatic Green function in vacuo:
Definition: pcm.F90:2146
pcm_oct_m::d_i_matrix
subroutine d_i_matrix(n_tess, tess)
Definition: pcm.F90:2126
pcm_oct_m::pcm_td_eom
integer, parameter, public pcm_td_eom
Definition: pcm.F90:275
pcm_oct_m::pcm_charges
subroutine, public pcm_charges(q_pcm, q_pcm_tot, v_cav, pcm_mat, n_tess, qtot_nominal, epsilon, renorm_charges, q_tot_tol, deltaQ)
Calculates the polarization charges at each tessera by using the response matrix 'pcm_mat',...
Definition: pcm.F90:1644
pcm_oct_m::pcm_vdw_species
integer, parameter pcm_vdw_species
Definition: pcm.F90:284
pcm_oct_m::pcm_v_cav_li
subroutine, public pcm_v_cav_li(v_cav, v_mesh, pcm, mesh)
Calculates the Hartree/external/kick potential at the tessera representative points by doing a 3D lin...
Definition: pcm.F90:1486
pcm_oct_m::pcm_calc_poisson
integer, parameter, public pcm_calc_poisson
Definition: pcm.F90:280
pcm_oct_m::gaubon
subroutine gaubon(sfe, nv, ns, pts, ccc, pp, pp1, area, intsph)
Use the Gauss-Bonnet theorem to calculate the area of the tessera with vertices 'pts(3,...
Definition: pcm.F90:2920
pcm_oct_m::mat_gamess
real(real64), dimension(:,:), allocatable mat_gamess
PCM matrix formatted to be inputed to GAMESS.
Definition: pcm.F90:273
pcm_oct_m::pcm_poisson_sanity_check
subroutine pcm_poisson_sanity_check(pcm, mesh)
Check that all the required nearest neighbors are prensent in the mesh.
Definition: pcm.F90:1729
pcm_oct_m::pcm_pot_rs_poisson
subroutine pcm_pot_rs_poisson(namespace, v_pcm, psolver, rho)
Generates the potential 'v_pcm' in real-space solving the poisson equation for rho.
Definition: pcm.F90:1926
pcm_oct_m::d_mat_elem_i
real(real64) function d_mat_elem_i(tessi, tessj)
Gradient of the Green function in vacuo .
Definition: pcm.F90:2178
pcm_oct_m::s_i_matrix
subroutine s_i_matrix(n_tess, tess)
Definition: pcm.F90:2107
pcm_oct_m::pcm_pot_rs_direct
subroutine pcm_pot_rs_direct(v_pcm, q_pcm, tess, n_tess, mesh, width_factor)
Generates the potential 'v_pcm' in real-space by direct sum.
Definition: pcm.F90:1942
pcm_oct_m::pcm_eps
subroutine, public pcm_eps(pcm, eps, omega)
Definition: pcm.F90:3255
pcm_oct_m::pcm_v_nuclei_cav
subroutine, public pcm_v_nuclei_cav(v_n_cav, ions, tess, n_tess)
Calculates the classical electrostatic potential geneated by the nuclei at the tesserae....
Definition: pcm.F90:1515
pcm_oct_m::pcm_end
subroutine, public pcm_end(pcm)
Definition: pcm.F90:3066
pcm_oct_m::pcm_vdw_optimized
integer, parameter pcm_vdw_optimized
Definition: pcm.F90:284
pcm_oct_m::pcm_init
subroutine, public pcm_init(pcm, namespace, space, ions, grid, qtot, val_charge, external_potentials_present, kick_present)
Initializes the PCM calculation: reads the VdW molecular cavity and generates the PCM response matrix...
Definition: pcm.F90:296
pcm_oct_m::pcm_td_neq
integer, parameter, public pcm_td_neq
Definition: pcm.F90:275
pcm_oct_m::vecp
subroutine vecp(p1, p2, p3, dnorm)
calculates the vectorial product p3 = p1 x p2
Definition: pcm.F90:3051
pcm_oct_m::inter
subroutine inter(sfe, p1, p2, p3, p4, ns, ia)
Finds the point 'p4', on the arc 'p1'-'p2' developed from 'p3', which is on the surface of sphere 'ns...
Definition: pcm.F90:2843
pcm_oct_m::pcm_td_eq
integer, parameter, public pcm_td_eq
Definition: pcm.F90:275
pcm_oct_m::subtessera
subroutine subtessera(sfe, ns, nesf, nv, pts, ccc, pp, pp1, area)
find the uncovered region for each tessera and computes the area, the representative point (pp) and t...
Definition: pcm.F90:2599
pcm_oct_m::pcm_min_input_parsing_for_spectrum
subroutine, public pcm_min_input_parsing_for_spectrum(pcm, namespace)
Definition: pcm.F90:3279
pcm_oct_m::cav_gen
subroutine cav_gen(tess_sphere, tess_min_distance, nesf, sfe, nts, cts, unit_pcminfo)
It builds the solute cavity surface and calculates the vertices, representative points and areas of t...
Definition: pcm.F90:2222
poisson_oct_m
Definition: poisson.F90:117
poisson_oct_m::dpoisson_solve
subroutine, public dpoisson_solve(this, namespace, pot, rho, all_nodes, kernel, reset)
Calculates the Poisson equation. Given the density returns the corresponding potential.
Definition: poisson.F90:867
profiling_oct_m
Definition: profiling.F90:118
profiling_oct_m::profiling_out
subroutine, public profiling_out(label)
Increment out counter and sum up difference between entry and exit time.
Definition: profiling.F90:631
profiling_oct_m::profiling_in
subroutine, public profiling_in(label, exclude)
Increment in counter and save entry time.
Definition: profiling.F90:554
space_oct_m
Definition: space.F90:116
species_oct_m
Definition: species.F90:116
unit_oct_m
brief This module defines the class unit_t which is used by the unit_systems_oct_m module.
Definition: unit.F90:134
unit_system_oct_m
This module defines the unit system, used for input and output.
Definition: unit_system.F90:128
unit_system_oct_m::units_out
type(unit_system_t), public units_out
Definition: unit_system.F90:163
unit_system_oct_m::units_inp
type(unit_system_t), public units_inp
the units systems for reading and writing
Definition: unit_system.F90:163
unit_system_oct_m::unit_one
type(unit_t), public unit_one
some special units required for particular quantities
Definition: unit_system.F90:166
varinfo_oct_m
Definition: varinfo.F90:116
box_minimum_oct_m::box_minimum_t
Class implementing a box that is a union of spheres. We do this in a specific class instead of using ...
Definition: box_minimum.F90:138
grid_oct_m::grid_t
Description of the grid, containing information on derivatives, stencil, and symmetries.
Definition: grid.F90:171
ions_oct_m::ions_t
Definition: ions.F90:164
mesh_interpolation_oct_m::mesh_interpolation_t
Definition: mesh_interpolation.F90:141
mesh_oct_m::mesh_t
Describes mesh distribution to nodes.
Definition: mesh.F90:187
namespace_oct_m::namespace_t
Definition: namespace.F90:117
pcm_eom_oct_m::pcm_tessera_t
tesselation derived type
Definition: pcm_eom.F90:144
pcm_oct_m::pcm_min_t
Definition: pcm.F90:257
pcm_oct_m::pcm_sphere_t
The cavity hosting the solute molecule is built from a set of interlocking spheres with optimized rad...
Definition: pcm.F90:176
pcm_oct_m::pcm_t
Definition: pcm.F90:187
poisson_oct_m::poisson_t
Definition: poisson.F90:199
space_oct_m::space_t
Definition: space.F90:132
true
int true(void)
Definition: symmetries_finite.c:3150