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