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