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