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