1!--------------------------------------------------------------------------------------------------!
2!   CP2K: A general program to perform molecular dynamics simulations                              !
3!   Copyright (C) 2000 - 2020  CP2K developers group                                               !
4!--------------------------------------------------------------------------------------------------!
5
6! **************************************************************************************************
7!> \brief Definition of the atomic potential types.
8!> \par History
9!>      GT, 22.09.2002: added elp_potential_types
10!> \author Matthias Krack (04.07.2000)
11! **************************************************************************************************
12MODULE external_potential_types
13   USE ao_util,                         ONLY: exp_radius
14   USE bibliography,                    ONLY: Goedecker1996,&
15                                              Hartwigsen1998,&
16                                              Krack2000,&
17                                              Krack2005,&
18                                              cite_reference
19   USE cp_linked_list_input,            ONLY: cp_sll_val_next,&
20                                              cp_sll_val_type
21   USE cp_para_types,                   ONLY: cp_para_env_type
22   USE cp_parser_methods,               ONLY: parser_get_next_line,&
23                                              parser_get_object,&
24                                              parser_search_string,&
25                                              parser_test_next_token
26   USE cp_parser_types,                 ONLY: cp_parser_type,&
27                                              parser_create,&
28                                              parser_release
29   USE input_section_types,             ONLY: section_vals_get,&
30                                              section_vals_list_get,&
31                                              section_vals_type,&
32                                              section_vals_val_set
33   USE input_val_types,                 ONLY: val_get,&
34                                              val_type
35   USE kinds,                           ONLY: default_path_length,&
36                                              default_string_length,&
37                                              dp
38   USE mathconstants,                   ONLY: dfac,&
39                                              fac,&
40                                              pi,&
41                                              rootpi
42   USE mathlib,                         ONLY: symmetrize_matrix
43   USE memory_utilities,                ONLY: reallocate
44   USE orbital_pointers,                ONLY: co,&
45                                              coset,&
46                                              init_orbital_pointers,&
47                                              nco,&
48                                              ncoset,&
49                                              nso
50   USE orbital_transformation_matrices, ONLY: orbtramat
51   USE periodic_table,                  ONLY: ptable
52   USE string_utilities,                ONLY: remove_word,&
53                                              uppercase
54#include "../base/base_uses.f90"
55
56   IMPLICIT NONE
57
58   PRIVATE
59
60   ! Global parameters
61
62   CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'external_potential_types'
63
64   ! Define the all-electron potential type
65
66   ! Literature: M. Krack and M. Parrinello,
67   !             Phys. Chem. Chem. Phys. 2, 2105 (2000)
68
69! **************************************************************************************************
70   TYPE all_potential_type
71      !MK PRIVATE
72      CHARACTER(LEN=default_string_length)   :: name
73      CHARACTER(LEN=default_string_length), &
74         DIMENSION(2)                      :: description
75      REAL(KIND=dp)                        :: alpha_core_charge, &
76                                              ccore_charge, &
77                                              core_charge_radius, &
78                                              zeff, zeff_correction
79      INTEGER                                :: z
80      INTEGER, DIMENSION(:), POINTER         :: elec_conf
81   END TYPE all_potential_type
82
83   ! Define the effective charge & inducible dipole potential type (for Fist)
84! **************************************************************************************************
85   TYPE fist_potential_type
86      PRIVATE
87      CHARACTER(LEN=default_string_length)     :: name
88      CHARACTER(LEN=default_string_length), &
89         DIMENSION(1)                        :: description
90      REAL(KIND=dp)                          :: apol, cpol, mm_radius, qeff, &
91                                                qmmm_corr_radius, qmmm_radius
92
93   END TYPE fist_potential_type
94
95! **************************************************************************************************
96! local potential type
97! V(r) = SUM_i exp(0.5*(r/rci)**2) * ( C1i + C2i (r/rci)**2 + C3i (r/rci)**4 ...)
98! alpha = 0.5/rci**2
99   TYPE local_potential_type
100      !PRIVATE
101      CHARACTER(LEN=default_string_length)       :: name
102      CHARACTER(LEN=default_string_length), &
103         DIMENSION(4)                          :: description
104      INTEGER                                    :: ngau, npol
105      REAL(KIND=dp)                            :: radius
106      REAL(KIND=dp), DIMENSION(:), POINTER     :: alpha
107      REAL(KIND=dp), DIMENSION(:, :), POINTER   :: cval
108   END TYPE local_potential_type
109
110   ! Define the GTH potential type
111
112   ! Literature: - S. Goedecker, M. Teter and J. Hutter,
113   !               Phys. Rev. B 54, 1703 (1996)
114   !             - C. Hartwigsen, S. Goedecker and J. Hutter,
115   !               Phys. Rev. B 58, 3641 (1998)
116   !             - M. Krack,
117   !               Theor. Chem. Acc. 114, 145 (2005)
118
119! **************************************************************************************************
120   TYPE gth_potential_type
121      !PRIVATE
122      CHARACTER(LEN=default_string_length)       :: name
123      CHARACTER(LEN=default_string_length)       :: aliases = ""
124      CHARACTER(LEN=default_string_length), &
125         DIMENSION(4)                          :: description
126      REAL(KIND=dp)                            :: alpha_core_charge, &
127                                                  alpha_ppl, ccore_charge, &
128                                                  cerf_ppl, zeff, &
129                                                  core_charge_radius, &
130                                                  ppl_radius, ppnl_radius, &
131                                                  zeff_correction
132      INTEGER                                    :: lppnl, lprj_ppnl_max, &
133                                                    nexp_ppl, nppnl, &
134                                                    nprj_ppnl_max, z
135      REAL(KIND=dp), DIMENSION(:), POINTER     :: alpha_ppnl, cexp_ppl
136      INTEGER, DIMENSION(:), POINTER             :: elec_conf
137      ! nonlocal projectors
138      INTEGER, DIMENSION(:), POINTER             :: nprj_ppnl
139      REAL(KIND=dp), DIMENSION(:, :), POINTER   :: cprj, cprj_ppnl, vprj_ppnl
140      REAL(KIND=dp), DIMENSION(:, :, :), POINTER :: hprj_ppnl
141      ! type extensions
142      ! NLCC
143      LOGICAL                                    :: nlcc
144      INTEGER                                    :: nexp_nlcc
145      REAL(KIND=dp), DIMENSION(:), POINTER     :: alpha_nlcc
146      INTEGER, DIMENSION(:), POINTER             :: nct_nlcc
147      REAL(KIND=dp), DIMENSION(:, :), POINTER   :: cval_nlcc
148      ! LSD potential
149      LOGICAL                                    :: lsdpot
150      INTEGER                                    :: nexp_lsd
151      REAL(KIND=dp), DIMENSION(:), POINTER     :: alpha_lsd
152      INTEGER, DIMENSION(:), POINTER             :: nct_lsd
153      REAL(KIND=dp), DIMENSION(:, :), POINTER   :: cval_lsd
154      ! extended local potential
155      LOGICAL                                    :: lpotextended
156      INTEGER                                    :: nexp_lpot
157      REAL(KIND=dp), DIMENSION(:), POINTER     :: alpha_lpot
158      INTEGER, DIMENSION(:), POINTER             :: nct_lpot
159      REAL(KIND=dp), DIMENSION(:, :), POINTER   :: cval_lpot
160
161   END TYPE gth_potential_type
162
163! **************************************************************************************************
164   TYPE sgp_potential_type
165      CHARACTER(LEN=default_string_length)       :: name
166      CHARACTER(LEN=default_string_length)       :: aliases = ""
167      CHARACTER(LEN=default_string_length), &
168         DIMENSION(4)                            :: description
169      ! CHARGE
170      INTEGER                                    :: z
171      REAL(KIND=dp)                              :: zeff, &
172                                                    zeff_correction
173      REAL(KIND=dp)                              :: alpha_core_charge, &
174                                                    ccore_charge, &
175                                                    core_charge_radius
176      REAL(KIND=dp)                              :: ppl_radius, ppnl_radius
177      INTEGER, DIMENSION(:), POINTER             :: elec_conf
178      ! LOCAL
179      LOGICAL                                    :: ecp_local
180      INTEGER                                    :: n_local
181      REAL(KIND=dp), DIMENSION(:), POINTER       :: a_local => Null()
182      REAL(KIND=dp), DIMENSION(:), POINTER       :: c_local => Null()
183      ! ECP local
184      INTEGER                                    :: nloc ! # terms
185      INTEGER, DIMENSION(1:10)                   :: nrloc ! r**(n-2)
186      REAL(dp), DIMENSION(1:10)                  :: aloc ! coefficient
187      REAL(dp), DIMENSION(1:10)                  :: bloc ! exponent
188      ! NONLOCAL
189      INTEGER                                    :: n_nonlocal
190      INTEGER                                    :: nppnl
191      INTEGER                                    :: lmax
192      LOGICAL, DIMENSION(0:5)                    :: is_nonlocal = .FALSE.
193      REAL(KIND=dp), DIMENSION(:), POINTER       :: a_nonlocal => Null()
194      REAL(KIND=dp), DIMENSION(:, :), POINTER    :: h_nonlocal => Null()
195      REAL(KIND=dp), DIMENSION(:, :, :), POINTER :: c_nonlocal => Null()
196      REAL(KIND=dp), DIMENSION(:, :), POINTER    :: cprj_ppnl
197      REAL(KIND=dp), DIMENSION(:), POINTER       :: vprj_ppnl
198      ! NLCC
199      LOGICAL                                    :: has_nlcc
200      INTEGER                                    :: n_nlcc
201      REAL(KIND=dp), DIMENSION(:), POINTER       :: a_nlcc => Null()
202      REAL(KIND=dp), DIMENSION(:), POINTER       :: c_nlcc => Null()
203   END TYPE sgp_potential_type
204
205! **************************************************************************************************
206   TYPE all_potential_p_type
207      TYPE(all_potential_type), POINTER          :: all_potential
208   END TYPE all_potential_p_type
209
210   TYPE gth_potential_p_type
211      TYPE(gth_potential_type), POINTER          :: gth_potential
212   END TYPE gth_potential_p_type
213
214   TYPE local_potential_p_type
215      TYPE(local_potential_type), POINTER        :: local_potential
216   END TYPE local_potential_p_type
217
218   TYPE sgp_potential_p_type
219      TYPE(sgp_potential_type), POINTER          :: sgp_potential
220   END TYPE sgp_potential_p_type
221! **************************************************************************************************
222
223   ! Public subroutines
224   PUBLIC :: allocate_potential, &
225             deallocate_potential, &
226             get_potential, &
227             init_potential, &
228             read_potential, &
229             set_potential, &
230             set_default_all_potential, &
231             write_potential, &
232             copy_potential
233
234   ! Public data types
235
236   PUBLIC :: all_potential_type, &
237             fist_potential_type, &
238             local_potential_type, &
239             gth_potential_type, &
240             sgp_potential_type
241   PUBLIC :: gth_potential_p_type, &
242             sgp_potential_p_type
243
244   INTERFACE allocate_potential
245      MODULE PROCEDURE allocate_all_potential, &
246         allocate_fist_potential, &
247         allocate_local_potential, &
248         allocate_gth_potential, &
249         allocate_sgp_potential
250   END INTERFACE
251
252   INTERFACE deallocate_potential
253      MODULE PROCEDURE deallocate_all_potential, &
254         deallocate_fist_potential, &
255         deallocate_local_potential, &
256         deallocate_sgp_potential, &
257         deallocate_gth_potential
258   END INTERFACE
259
260   INTERFACE get_potential
261      MODULE PROCEDURE get_all_potential, &
262         get_fist_potential, &
263         get_local_potential, &
264         get_gth_potential, &
265         get_sgp_potential
266   END INTERFACE
267
268   INTERFACE init_potential
269      MODULE PROCEDURE init_all_potential, &
270         init_gth_potential, &
271         init_sgp_potential
272   END INTERFACE
273
274   INTERFACE read_potential
275      MODULE PROCEDURE read_all_potential, &
276         read_local_potential, &
277         read_gth_potential
278   END INTERFACE
279
280   INTERFACE set_potential
281      MODULE PROCEDURE set_all_potential, &
282         set_fist_potential, &
283         set_local_potential, &
284         set_gth_potential, &
285         set_sgp_potential
286   END INTERFACE
287
288   INTERFACE write_potential
289      MODULE PROCEDURE write_all_potential, &
290         write_local_potential, &
291         write_gth_potential, &
292         write_sgp_potential
293   END INTERFACE
294
295   INTERFACE copy_potential
296      MODULE PROCEDURE copy_all_potential, &
297         copy_gth_potential, &
298         copy_sgp_potential
299   END INTERFACE
300
301CONTAINS
302
303! **************************************************************************************************
304!> \brief   Allocate an atomic all-electron potential data set.
305!> \param potential ...
306!> \date    25.07.2000,
307!> \author  MK
308!> \version 1.0
309! **************************************************************************************************
310   SUBROUTINE allocate_all_potential(potential)
311      TYPE(all_potential_type), POINTER                  :: potential
312
313      CHARACTER(len=*), PARAMETER :: routineN = 'allocate_all_potential', &
314         routineP = moduleN//':'//routineN
315
316      IF (ASSOCIATED(potential)) CALL deallocate_potential(potential)
317
318      ALLOCATE (potential)
319
320      NULLIFY (potential%elec_conf)
321
322      potential%description(1) = "All-electron potential"
323      potential%description(2) = "Krack, Parrinello, PCCP 2, 2105 (2000)"
324
325   END SUBROUTINE allocate_all_potential
326
327! **************************************************************************************************
328!> \brief   Allocate an effective charge and inducible dipole potential data set.
329!> \param potential ...
330!> \date    05.03.2010
331!> \author  Toon.Verstraelen@gmail.com
332! **************************************************************************************************
333   SUBROUTINE allocate_fist_potential(potential)
334      TYPE(fist_potential_type), POINTER                 :: potential
335
336      CHARACTER(len=*), PARAMETER :: routineN = 'allocate_fist_potential', &
337         routineP = moduleN//':'//routineN
338
339      IF (ASSOCIATED(potential)) CALL deallocate_potential(potential)
340
341      ALLOCATE (potential)
342
343      potential%apol = 0.0_dp
344      potential%cpol = 0.0_dp
345      potential%mm_radius = 0.0_dp
346      potential%qeff = 0.0_dp
347      potential%qmmm_radius = 0.0_dp
348      potential%qmmm_corr_radius = 0.0_dp
349
350      potential%description(1) = "Effective charge and inducible dipole potential"
351
352   END SUBROUTINE allocate_fist_potential
353
354! **************************************************************************************************
355!> \brief   Allocate an atomic local potential data set.
356!> \param potential ...
357!> \date    24.01.2014
358!> \author  JGH
359!> \version 1.0
360! **************************************************************************************************
361   SUBROUTINE allocate_local_potential(potential)
362      TYPE(local_potential_type), POINTER                :: potential
363
364      CHARACTER(len=*), PARAMETER :: routineN = 'allocate_local_potential', &
365         routineP = moduleN//':'//routineN
366
367      IF (ASSOCIATED(potential)) CALL deallocate_potential(potential)
368
369      ALLOCATE (potential)
370
371      NULLIFY (potential%alpha)
372      NULLIFY (potential%cval)
373
374      potential%description(1) = "Local short-range pseudopotential"
375      potential%ngau = 0
376      potential%npol = 0
377      potential%radius = 0.0_dp
378
379   END SUBROUTINE allocate_local_potential
380
381! **************************************************************************************************
382!> \brief   Allocate an atomic GTH potential data set.
383!> \param potential ...
384!> \date    25.07.2000
385!> \author  MK
386!> \version 1.0
387! **************************************************************************************************
388   SUBROUTINE allocate_gth_potential(potential)
389      TYPE(gth_potential_type), POINTER                  :: potential
390
391      CHARACTER(len=*), PARAMETER :: routineN = 'allocate_gth_potential', &
392         routineP = moduleN//':'//routineN
393
394      IF (ASSOCIATED(potential)) CALL deallocate_potential(potential)
395
396      ALLOCATE (potential)
397
398      NULLIFY (potential%alpha_ppnl)
399      NULLIFY (potential%cexp_ppl)
400      NULLIFY (potential%elec_conf)
401      NULLIFY (potential%nprj_ppnl)
402      NULLIFY (potential%cprj)
403      NULLIFY (potential%cprj_ppnl)
404      NULLIFY (potential%vprj_ppnl)
405      NULLIFY (potential%hprj_ppnl)
406
407      NULLIFY (potential%alpha_lpot)
408      NULLIFY (potential%nct_lpot)
409      NULLIFY (potential%cval_lpot)
410      NULLIFY (potential%alpha_lsd)
411      NULLIFY (potential%nct_lsd)
412      NULLIFY (potential%cval_lsd)
413      NULLIFY (potential%alpha_nlcc)
414      NULLIFY (potential%nct_nlcc)
415      NULLIFY (potential%cval_nlcc)
416
417      potential%description(1) = "Goedecker-Teter-Hutter pseudopotential"
418      potential%description(2) = "Goedecker et al., PRB 54, 1703 (1996)"
419      potential%description(3) = "Hartwigsen et al., PRB 58, 3641 (1998)"
420      potential%description(4) = "Krack, TCA 114, 145 (2005)"
421
422   END SUBROUTINE allocate_gth_potential
423
424! **************************************************************************************************
425!> \brief   Allocate an atomic SGP potential data set.
426!> \param potential ...
427!> \version 1.0
428! **************************************************************************************************
429   SUBROUTINE allocate_sgp_potential(potential)
430      TYPE(sgp_potential_type), POINTER                  :: potential
431
432      CHARACTER(len=*), PARAMETER :: routineN = 'allocate_sgp_potential', &
433         routineP = moduleN//':'//routineN
434
435      IF (ASSOCIATED(potential)) CALL deallocate_potential(potential)
436
437      ALLOCATE (potential)
438
439      NULLIFY (potential%elec_conf)
440      NULLIFY (potential%a_local)
441      NULLIFY (potential%c_local)
442      NULLIFY (potential%a_nonlocal)
443      NULLIFY (potential%c_nonlocal)
444      NULLIFY (potential%h_nonlocal)
445      NULLIFY (potential%a_nlcc)
446      NULLIFY (potential%c_nlcc)
447      NULLIFY (potential%cprj_ppnl)
448      NULLIFY (potential%vprj_ppnl)
449
450      potential%name = ""
451      potential%aliases = ""
452      potential%description(1) = "Separable Gaussian pseudopotential"
453      potential%description(2) = "M. Pelissier, N. Komiha, J.P. Daudey, JCC, 9, 298 (1988)"
454      potential%description(3) = "create from"
455      potential%description(4) = ""
456
457      potential%z = 0
458      potential%zeff = 0.0_dp
459      potential%zeff_correction = 0.0_dp
460      potential%alpha_core_charge = 0.0_dp
461      potential%ccore_charge = 0.0_dp
462      potential%core_charge_radius = 0.0_dp
463      potential%ppl_radius = 0.0_dp
464      potential%ppnl_radius = 0.0_dp
465      potential%ecp_local = .FALSE.
466      potential%n_local = 0
467      potential%nloc = 0
468      potential%nrloc = 0
469      potential%aloc = 0.0_dp
470      potential%bloc = 0.0_dp
471      potential%n_nonlocal = 0
472      potential%nppnl = 0
473      potential%lmax = -1
474      potential%is_nonlocal = .FALSE.
475      potential%has_nlcc = .FALSE.
476      potential%n_nlcc = 0
477
478   END SUBROUTINE allocate_sgp_potential
479! **************************************************************************************************
480!> \brief   Deallocate an atomic all-electron potential data set.
481!> \param potential ...
482!> \date    03.11.2000
483!> \author  MK
484!> \version 1.0
485! **************************************************************************************************
486   SUBROUTINE deallocate_all_potential(potential)
487      TYPE(all_potential_type), POINTER                  :: potential
488
489      CHARACTER(len=*), PARAMETER :: routineN = 'deallocate_all_potential', &
490         routineP = moduleN//':'//routineN
491
492      IF (ASSOCIATED(potential)) THEN
493         DEALLOCATE (potential%elec_conf)
494         DEALLOCATE (potential)
495      ELSE
496         CPABORT("The pointer potential is not associated.")
497      END IF
498
499   END SUBROUTINE deallocate_all_potential
500
501! **************************************************************************************************
502!> \brief   Deallocate an effective charge and inducible dipole potential data set.
503!> \param potential ...
504!> \date    05.03.2010
505!> \author  Toon.Verstraelen@gmail.com
506! **************************************************************************************************
507   SUBROUTINE deallocate_fist_potential(potential)
508      TYPE(fist_potential_type), POINTER                 :: potential
509
510      CHARACTER(len=*), PARAMETER :: routineN = 'deallocate_fist_potential', &
511         routineP = moduleN//':'//routineN
512
513      IF (ASSOCIATED(potential)) THEN
514         ! Nothing exciting here yet.
515         DEALLOCATE (potential)
516      ELSE
517         CPABORT("The pointer potential is not associated.")
518      END IF
519
520   END SUBROUTINE deallocate_fist_potential
521
522! **************************************************************************************************
523!> \brief   Deallocate an atomic local potential data set.
524!> \param potential ...
525!> \date    24.01.2014
526!> \author  JGH
527!> \version 1.0
528! **************************************************************************************************
529   SUBROUTINE deallocate_local_potential(potential)
530      TYPE(local_potential_type), POINTER                :: potential
531
532      CHARACTER(len=*), PARAMETER :: routineN = 'deallocate_local_potential', &
533         routineP = moduleN//':'//routineN
534
535      IF (ASSOCIATED(potential)) THEN
536
537         IF (ASSOCIATED(potential%alpha)) THEN
538            DEALLOCATE (potential%alpha)
539         END IF
540         IF (ASSOCIATED(potential%cval)) THEN
541            DEALLOCATE (potential%cval)
542         END IF
543
544         DEALLOCATE (potential)
545
546      ELSE
547
548         CPABORT("The pointer potential is not associated.")
549
550      END IF
551
552   END SUBROUTINE deallocate_local_potential
553
554! **************************************************************************************************
555!> \brief   Deallocate an atomic GTH potential data set.
556!> \param potential ...
557!> \date    03.11.2000
558!> \author  MK
559!> \version 1.0
560! **************************************************************************************************
561   SUBROUTINE deallocate_gth_potential(potential)
562      TYPE(gth_potential_type), POINTER                  :: potential
563
564      CHARACTER(len=*), PARAMETER :: routineN = 'deallocate_gth_potential', &
565         routineP = moduleN//':'//routineN
566
567      IF (ASSOCIATED(potential)) THEN
568
569         DEALLOCATE (potential%elec_conf)
570         !     *** Deallocate the parameters of the local part ***
571
572         IF (ASSOCIATED(potential%cexp_ppl)) THEN
573            DEALLOCATE (potential%cexp_ppl)
574         END IF
575
576         !     *** Deallocate the parameters of the non-local part ***
577         IF (ASSOCIATED(potential%alpha_ppnl)) THEN
578            DEALLOCATE (potential%alpha_ppnl)
579            DEALLOCATE (potential%cprj)
580            DEALLOCATE (potential%cprj_ppnl)
581            DEALLOCATE (potential%hprj_ppnl)
582            DEALLOCATE (potential%nprj_ppnl)
583            DEALLOCATE (potential%vprj_ppnl)
584         END IF
585
586         IF (ASSOCIATED(potential%alpha_lpot)) THEN
587            DEALLOCATE (potential%alpha_lpot)
588            DEALLOCATE (potential%nct_lpot)
589            DEALLOCATE (potential%cval_lpot)
590         END IF
591
592         IF (ASSOCIATED(potential%alpha_lsd)) THEN
593            DEALLOCATE (potential%alpha_lsd)
594            DEALLOCATE (potential%nct_lsd)
595            DEALLOCATE (potential%cval_lsd)
596         END IF
597
598         IF (ASSOCIATED(potential%alpha_nlcc)) THEN
599            DEALLOCATE (potential%alpha_nlcc)
600            DEALLOCATE (potential%nct_nlcc)
601            DEALLOCATE (potential%cval_nlcc)
602         END IF
603
604         DEALLOCATE (potential)
605      ELSE
606         CPABORT("The pointer potential is not associated.")
607      END IF
608
609   END SUBROUTINE deallocate_gth_potential
610
611! **************************************************************************************************
612!> \brief   Deallocate an atomic SGP potential data set.
613!> \param potential ...
614! **************************************************************************************************
615   SUBROUTINE deallocate_sgp_potential(potential)
616      TYPE(sgp_potential_type), POINTER                  :: potential
617
618      CHARACTER(len=*), PARAMETER :: routineN = 'deallocate_sgp_potential', &
619         routineP = moduleN//':'//routineN
620
621      IF (ASSOCIATED(potential)) THEN
622         IF (ASSOCIATED(potential%elec_conf)) THEN
623            DEALLOCATE (potential%elec_conf)
624         END IF
625         IF (ASSOCIATED(potential%a_local)) THEN
626            DEALLOCATE (potential%a_local)
627         END IF
628         IF (ASSOCIATED(potential%c_local)) THEN
629            DEALLOCATE (potential%c_local)
630         END IF
631
632         IF (ASSOCIATED(potential%a_nonlocal)) THEN
633            DEALLOCATE (potential%a_nonlocal)
634         END IF
635         IF (ASSOCIATED(potential%h_nonlocal)) THEN
636            DEALLOCATE (potential%h_nonlocal)
637         END IF
638         IF (ASSOCIATED(potential%c_nonlocal)) THEN
639            DEALLOCATE (potential%c_nonlocal)
640         END IF
641         IF (ASSOCIATED(potential%cprj_ppnl)) THEN
642            DEALLOCATE (potential%cprj_ppnl)
643         END IF
644         IF (ASSOCIATED(potential%vprj_ppnl)) THEN
645            DEALLOCATE (potential%vprj_ppnl)
646         END IF
647
648         IF (ASSOCIATED(potential%a_nlcc)) THEN
649            DEALLOCATE (potential%a_nlcc)
650         END IF
651         IF (ASSOCIATED(potential%c_nlcc)) THEN
652            DEALLOCATE (potential%c_nlcc)
653         END IF
654
655         DEALLOCATE (potential)
656      ELSE
657         CPABORT("The pointer potential is not associated.")
658      END IF
659
660   END SUBROUTINE deallocate_sgp_potential
661
662! **************************************************************************************************
663!> \brief   Get attributes of an all-electron potential data set.
664!> \param potential ...
665!> \param name ...
666!> \param alpha_core_charge ...
667!> \param ccore_charge ...
668!> \param core_charge_radius ...
669!> \param z ...
670!> \param zeff ...
671!> \param zeff_correction ...
672!> \param elec_conf ...
673!> \date    11.01.2002
674!> \author  MK
675!> \version 1.0
676! **************************************************************************************************
677   SUBROUTINE get_all_potential(potential, name, alpha_core_charge, &
678                                ccore_charge, core_charge_radius, z, zeff, &
679                                zeff_correction, elec_conf)
680      TYPE(all_potential_type), POINTER                  :: potential
681      CHARACTER(LEN=default_string_length), &
682         INTENT(OUT), OPTIONAL                           :: name
683      REAL(KIND=dp), INTENT(OUT), OPTIONAL               :: alpha_core_charge, ccore_charge, &
684                                                            core_charge_radius
685      INTEGER, INTENT(OUT), OPTIONAL                     :: z
686      REAL(KIND=dp), INTENT(OUT), OPTIONAL               :: zeff, zeff_correction
687      INTEGER, DIMENSION(:), OPTIONAL, POINTER           :: elec_conf
688
689      CHARACTER(len=*), PARAMETER :: routineN = 'get_all_potential', &
690         routineP = moduleN//':'//routineN
691
692      IF (ASSOCIATED(potential)) THEN
693
694         IF (PRESENT(name)) name = potential%name
695         IF (PRESENT(alpha_core_charge)) &
696            alpha_core_charge = potential%alpha_core_charge
697         IF (PRESENT(ccore_charge)) ccore_charge = potential%ccore_charge
698         IF (PRESENT(core_charge_radius)) &
699            core_charge_radius = potential%core_charge_radius
700         IF (PRESENT(z)) z = potential%z
701         IF (PRESENT(zeff)) zeff = potential%zeff
702         IF (PRESENT(zeff_correction)) zeff_correction = potential%zeff_correction
703         IF (PRESENT(elec_conf)) elec_conf => potential%elec_conf
704
705      ELSE
706
707         CPABORT("The pointer potential is not associated.")
708
709      END IF
710
711   END SUBROUTINE get_all_potential
712
713! **************************************************************************************************
714!> \brief   Get attributes of an effective point charge and inducible dipole
715!>          potential.
716!> \param potential ...
717!> \param name ...
718!> \param apol ...
719!> \param cpol ...
720!> \param mm_radius ...
721!> \param qeff ...
722!> \param qmmm_corr_radius ...
723!> \param qmmm_radius ...
724!> \date    05.03-2010
725!> \author  Toon.Verstraelen@UGent.be
726! **************************************************************************************************
727   SUBROUTINE get_fist_potential(potential, name, apol, cpol, mm_radius, qeff, &
728                                 qmmm_corr_radius, qmmm_radius)
729      TYPE(fist_potential_type), POINTER                 :: potential
730      CHARACTER(LEN=default_string_length), &
731         INTENT(OUT), OPTIONAL                           :: name
732      REAL(KIND=dp), INTENT(OUT), OPTIONAL               :: apol, cpol, mm_radius, qeff, &
733                                                            qmmm_corr_radius, qmmm_radius
734
735      CHARACTER(len=*), PARAMETER :: routineN = 'get_fist_potential', &
736         routineP = moduleN//':'//routineN
737
738      IF (ASSOCIATED(potential)) THEN
739
740         IF (PRESENT(name)) name = potential%name
741         IF (PRESENT(apol)) apol = potential%apol
742         IF (PRESENT(cpol)) cpol = potential%cpol
743         IF (PRESENT(mm_radius)) mm_radius = potential%mm_radius
744         IF (PRESENT(qeff)) qeff = potential%qeff
745         IF (PRESENT(qmmm_corr_radius)) qmmm_corr_radius = potential%qmmm_corr_radius
746         IF (PRESENT(qmmm_radius)) qmmm_radius = potential%qmmm_radius
747
748      ELSE
749
750         CPABORT("The pointer potential is not associated.")
751
752      END IF
753
754   END SUBROUTINE get_fist_potential
755
756! **************************************************************************************************
757!> \brief   Get attributes of an atomic local potential data set.
758!> \param potential ...
759!> \param name ...
760!> \param ngau ...
761!> \param npol ...
762!> \param alpha ...
763!> \param cval ...
764!> \param radius ...
765!> \date    24.01.2014
766!> \author  JGH
767!> \version 1.0
768! **************************************************************************************************
769   SUBROUTINE get_local_potential(potential, name, ngau, npol, alpha, cval, radius)
770      TYPE(local_potential_type), POINTER                :: potential
771      CHARACTER(LEN=default_string_length), &
772         INTENT(OUT), OPTIONAL                           :: name
773      INTEGER, INTENT(OUT), OPTIONAL                     :: ngau, npol
774      REAL(KIND=dp), DIMENSION(:), OPTIONAL, POINTER     :: alpha
775      REAL(KIND=dp), DIMENSION(:, :), OPTIONAL, POINTER  :: cval
776      REAL(KIND=dp), INTENT(OUT), OPTIONAL               :: radius
777
778      CHARACTER(len=*), PARAMETER :: routineN = 'get_local_potential', &
779         routineP = moduleN//':'//routineN
780
781      IF (ASSOCIATED(potential)) THEN
782
783         IF (PRESENT(name)) name = potential%name
784         IF (PRESENT(ngau)) ngau = potential%ngau
785         IF (PRESENT(npol)) npol = potential%npol
786         IF (PRESENT(alpha)) alpha => potential%alpha
787         IF (PRESENT(cval)) cval => potential%cval
788         IF (PRESENT(radius)) radius = potential%radius
789
790      ELSE
791
792         CPABORT("The pointer potential is not associated.")
793
794      END IF
795
796   END SUBROUTINE get_local_potential
797
798! **************************************************************************************************
799!> \brief   Get attributes of a GTH potential data set.
800!> \param potential ...
801!> \param name ...
802!> \param aliases ...
803!> \param alpha_core_charge ...
804!> \param alpha_ppl ...
805!> \param ccore_charge ...
806!> \param cerf_ppl ...
807!> \param core_charge_radius ...
808!> \param ppl_radius ...
809!> \param ppnl_radius ...
810!> \param lppnl ...
811!> \param lprj_ppnl_max ...
812!> \param nexp_ppl ...
813!> \param nppnl ...
814!> \param nprj_ppnl_max ...
815!> \param z ...
816!> \param zeff ...
817!> \param zeff_correction ...
818!> \param ppl_present ...
819!> \param ppnl_present ...
820!> \param alpha_ppnl ...
821!> \param cexp_ppl ...
822!> \param elec_conf ...
823!> \param nprj_ppnl ...
824!> \param cprj ...
825!> \param cprj_ppnl ...
826!> \param vprj_ppnl ...
827!> \param hprj_ppnl ...
828!> \param lpot_present ...
829!> \param nexp_lpot ...
830!> \param alpha_lpot ...
831!> \param nct_lpot ...
832!> \param cval_lpot ...
833!> \param lsd_present ...
834!> \param nexp_lsd ...
835!> \param alpha_lsd ...
836!> \param nct_lsd ...
837!> \param cval_lsd ...
838!> \param nlcc_present ...
839!> \param nexp_nlcc ...
840!> \param alpha_nlcc ...
841!> \param nct_nlcc ...
842!> \param cval_nlcc ...
843!> \date    11.01.2002
844!> \author  MK
845!> \version 1.0
846! **************************************************************************************************
847   SUBROUTINE get_gth_potential(potential, name, aliases, alpha_core_charge, &
848                                alpha_ppl, ccore_charge, cerf_ppl, &
849                                core_charge_radius, ppl_radius, ppnl_radius, &
850                                lppnl, lprj_ppnl_max, nexp_ppl, nppnl, &
851                                nprj_ppnl_max, z, zeff, zeff_correction, &
852                                ppl_present, ppnl_present, &
853                                alpha_ppnl, cexp_ppl, elec_conf, nprj_ppnl, cprj, &
854                                cprj_ppnl, vprj_ppnl, hprj_ppnl, &
855                                lpot_present, nexp_lpot, alpha_lpot, nct_lpot, cval_lpot, &
856                                lsd_present, nexp_lsd, alpha_lsd, nct_lsd, cval_lsd, &
857                                nlcc_present, nexp_nlcc, alpha_nlcc, nct_nlcc, cval_nlcc)
858
859      TYPE(gth_potential_type), POINTER                  :: potential
860      CHARACTER(LEN=default_string_length), &
861         INTENT(OUT), OPTIONAL                           :: name, aliases
862      REAL(KIND=dp), INTENT(OUT), OPTIONAL               :: alpha_core_charge, alpha_ppl, &
863                                                            ccore_charge, cerf_ppl, &
864                                                            core_charge_radius, ppl_radius, &
865                                                            ppnl_radius
866      INTEGER, INTENT(OUT), OPTIONAL                     :: lppnl, lprj_ppnl_max, nexp_ppl, nppnl, &
867                                                            nprj_ppnl_max, z
868      REAL(KIND=dp), INTENT(OUT), OPTIONAL               :: zeff, zeff_correction
869      LOGICAL, INTENT(OUT), OPTIONAL                     :: ppl_present, ppnl_present
870      REAL(KIND=dp), DIMENSION(:), OPTIONAL, POINTER     :: alpha_ppnl, cexp_ppl
871      INTEGER, DIMENSION(:), OPTIONAL, POINTER           :: elec_conf, nprj_ppnl
872      REAL(KIND=dp), DIMENSION(:, :), OPTIONAL, POINTER  :: cprj, cprj_ppnl, vprj_ppnl
873      REAL(KIND=dp), DIMENSION(:, :, :), OPTIONAL, &
874         POINTER                                         :: hprj_ppnl
875      LOGICAL, INTENT(OUT), OPTIONAL                     :: lpot_present
876      INTEGER, INTENT(OUT), OPTIONAL                     :: nexp_lpot
877      REAL(KIND=dp), DIMENSION(:), OPTIONAL, POINTER     :: alpha_lpot
878      INTEGER, DIMENSION(:), OPTIONAL, POINTER           :: nct_lpot
879      REAL(KIND=dp), DIMENSION(:, :), OPTIONAL, POINTER  :: cval_lpot
880      LOGICAL, INTENT(OUT), OPTIONAL                     :: lsd_present
881      INTEGER, INTENT(OUT), OPTIONAL                     :: nexp_lsd
882      REAL(KIND=dp), DIMENSION(:), OPTIONAL, POINTER     :: alpha_lsd
883      INTEGER, DIMENSION(:), OPTIONAL, POINTER           :: nct_lsd
884      REAL(KIND=dp), DIMENSION(:, :), OPTIONAL, POINTER  :: cval_lsd
885      LOGICAL, INTENT(OUT), OPTIONAL                     :: nlcc_present
886      INTEGER, INTENT(OUT), OPTIONAL                     :: nexp_nlcc
887      REAL(KIND=dp), DIMENSION(:), OPTIONAL, POINTER     :: alpha_nlcc
888      INTEGER, DIMENSION(:), OPTIONAL, POINTER           :: nct_nlcc
889      REAL(KIND=dp), DIMENSION(:, :), OPTIONAL, POINTER  :: cval_nlcc
890
891      CHARACTER(len=*), PARAMETER :: routineN = 'get_gth_potential', &
892         routineP = moduleN//':'//routineN
893
894      IF (ASSOCIATED(potential)) THEN
895
896         IF (PRESENT(name)) name = potential%name
897         IF (PRESENT(aliases)) aliases = potential%aliases
898         IF (PRESENT(alpha_core_charge)) &
899            alpha_core_charge = potential%alpha_core_charge
900         IF (PRESENT(alpha_ppl)) alpha_ppl = potential%alpha_ppl
901         IF (PRESENT(ccore_charge)) ccore_charge = potential%ccore_charge
902         IF (PRESENT(cerf_ppl)) cerf_ppl = potential%cerf_ppl
903         IF (PRESENT(core_charge_radius)) &
904            core_charge_radius = potential%core_charge_radius
905         IF (PRESENT(ppl_radius)) ppl_radius = potential%ppl_radius
906         IF (PRESENT(ppnl_radius)) ppnl_radius = potential%ppnl_radius
907         IF (PRESENT(lppnl)) lppnl = potential%lppnl
908         IF (PRESENT(lprj_ppnl_max)) lprj_ppnl_max = potential%lprj_ppnl_max
909         IF (PRESENT(nexp_ppl)) nexp_ppl = potential%nexp_ppl
910         IF (PRESENT(nppnl)) nppnl = potential%nppnl
911         IF (PRESENT(nprj_ppnl_max)) nprj_ppnl_max = potential%nprj_ppnl_max
912         IF (PRESENT(z)) z = potential%z
913         IF (PRESENT(zeff)) zeff = potential%zeff
914         IF (PRESENT(zeff_correction)) zeff_correction = potential%zeff_correction
915         IF (PRESENT(ppl_present)) ppl_present = (potential%nexp_ppl > 0)
916         IF (PRESENT(ppnl_present)) ppnl_present = (potential%nppnl > 0)
917         IF (PRESENT(alpha_ppnl)) alpha_ppnl => potential%alpha_ppnl
918         IF (PRESENT(cexp_ppl)) cexp_ppl => potential%cexp_ppl
919         IF (PRESENT(elec_conf)) elec_conf => potential%elec_conf
920         IF (PRESENT(nprj_ppnl)) nprj_ppnl => potential%nprj_ppnl
921         IF (PRESENT(cprj)) cprj => potential%cprj
922         IF (PRESENT(cprj_ppnl)) cprj_ppnl => potential%cprj_ppnl
923         IF (PRESENT(vprj_ppnl)) vprj_ppnl => potential%vprj_ppnl
924         IF (PRESENT(hprj_ppnl)) hprj_ppnl => potential%hprj_ppnl
925
926         IF (PRESENT(lpot_present)) lpot_present = potential%lpotextended
927         IF (PRESENT(nexp_lpot)) nexp_lpot = potential%nexp_lpot
928         IF (PRESENT(alpha_lpot)) alpha_lpot => potential%alpha_lpot
929         IF (PRESENT(nct_lpot)) nct_lpot => potential%nct_lpot
930         IF (PRESENT(cval_lpot)) cval_lpot => potential%cval_lpot
931
932         IF (PRESENT(lsd_present)) lsd_present = potential%lsdpot
933         IF (PRESENT(nexp_lsd)) nexp_lsd = potential%nexp_lsd
934         IF (PRESENT(alpha_lsd)) alpha_lsd => potential%alpha_lsd
935         IF (PRESENT(nct_lsd)) nct_lsd => potential%nct_lsd
936         IF (PRESENT(cval_lsd)) cval_lsd => potential%cval_lsd
937
938         IF (PRESENT(nlcc_present)) nlcc_present = potential%nlcc
939         IF (PRESENT(nexp_nlcc)) nexp_nlcc = potential%nexp_nlcc
940         IF (PRESENT(alpha_nlcc)) alpha_nlcc => potential%alpha_nlcc
941         IF (PRESENT(nct_nlcc)) nct_nlcc => potential%nct_nlcc
942         IF (PRESENT(cval_nlcc)) cval_nlcc => potential%cval_nlcc
943
944      ELSE
945
946         CPABORT("The pointer potential is not associated.")
947
948      END IF
949
950   END SUBROUTINE get_gth_potential
951
952! **************************************************************************************************
953!> \brief ...
954!> \param potential ...
955!> \param name ...
956!> \param description ...
957!> \param aliases ...
958!> \param elec_conf ...
959!> \param z ...
960!> \param zeff ...
961!> \param zeff_correction ...
962!> \param alpha_core_charge ...
963!> \param ccore_charge ...
964!> \param core_charge_radius ...
965!> \param ppl_radius ...
966!> \param ppnl_radius ...
967!> \param ppl_present ...
968!> \param ppnl_present ...
969!> \param ecp_local ...
970!> \param n_local ...
971!> \param a_local ...
972!> \param c_local ...
973!> \param nloc ...
974!> \param nrloc ...
975!> \param aloc ...
976!> \param bloc ...
977!> \param n_nonlocal ...
978!> \param nppnl ...
979!> \param lmax ...
980!> \param is_nonlocal ...
981!> \param a_nonlocal ...
982!> \param h_nonlocal ...
983!> \param c_nonlocal ...
984!> \param cprj_ppnl ...
985!> \param vprj_ppnl ...
986!> \param has_nlcc ...
987!> \param n_nlcc ...
988!> \param a_nlcc ...
989!> \param c_nlcc ...
990! **************************************************************************************************
991   SUBROUTINE get_sgp_potential(potential, name, description, aliases, elec_conf, &
992                                z, zeff, zeff_correction, alpha_core_charge, &
993                                ccore_charge, core_charge_radius, &
994                                ppl_radius, ppnl_radius, ppl_present, ppnl_present, &
995                                ecp_local, n_local, a_local, c_local, &
996                                nloc, nrloc, aloc, bloc, &
997                                n_nonlocal, nppnl, lmax, is_nonlocal, a_nonlocal, h_nonlocal, c_nonlocal, &
998                                cprj_ppnl, vprj_ppnl, has_nlcc, n_nlcc, a_nlcc, c_nlcc)
999
1000      TYPE(sgp_potential_type), POINTER                  :: potential
1001      CHARACTER(LEN=default_string_length), &
1002         INTENT(OUT), OPTIONAL                           :: name
1003      CHARACTER(LEN=default_string_length), &
1004         DIMENSION(4), INTENT(OUT), OPTIONAL             :: description
1005      CHARACTER(LEN=default_string_length), &
1006         INTENT(OUT), OPTIONAL                           :: aliases
1007      INTEGER, DIMENSION(:), OPTIONAL, POINTER           :: elec_conf
1008      INTEGER, INTENT(OUT), OPTIONAL                     :: z
1009      REAL(KIND=dp), INTENT(OUT), OPTIONAL               :: zeff, zeff_correction, &
1010                                                            alpha_core_charge, ccore_charge, &
1011                                                            core_charge_radius, ppl_radius, &
1012                                                            ppnl_radius
1013      LOGICAL, INTENT(OUT), OPTIONAL                     :: ppl_present, ppnl_present, ecp_local
1014      INTEGER, INTENT(OUT), OPTIONAL                     :: n_local
1015      REAL(KIND=dp), DIMENSION(:), OPTIONAL, POINTER     :: a_local, c_local
1016      INTEGER, INTENT(OUT), OPTIONAL                     :: nloc
1017      INTEGER, DIMENSION(1:10), INTENT(OUT), OPTIONAL    :: nrloc
1018      REAL(dp), DIMENSION(1:10), INTENT(OUT), OPTIONAL   :: aloc, bloc
1019      INTEGER, INTENT(OUT), OPTIONAL                     :: n_nonlocal, nppnl, lmax
1020      LOGICAL, DIMENSION(0:5), OPTIONAL                  :: is_nonlocal
1021      REAL(KIND=dp), DIMENSION(:), OPTIONAL, POINTER     :: a_nonlocal
1022      REAL(KIND=dp), DIMENSION(:, :), OPTIONAL, POINTER  :: h_nonlocal
1023      REAL(KIND=dp), DIMENSION(:, :, :), OPTIONAL, &
1024         POINTER                                         :: c_nonlocal
1025      REAL(KIND=dp), DIMENSION(:, :), OPTIONAL, POINTER  :: cprj_ppnl
1026      REAL(KIND=dp), DIMENSION(:), OPTIONAL, POINTER     :: vprj_ppnl
1027      LOGICAL, INTENT(OUT), OPTIONAL                     :: has_nlcc
1028      INTEGER, INTENT(OUT), OPTIONAL                     :: n_nlcc
1029      REAL(KIND=dp), DIMENSION(:), OPTIONAL, POINTER     :: a_nlcc, c_nlcc
1030
1031      CHARACTER(len=*), PARAMETER :: routineN = 'get_sgp_potential', &
1032         routineP = moduleN//':'//routineN
1033
1034      IF (ASSOCIATED(potential)) THEN
1035
1036         IF (PRESENT(name)) name = potential%name
1037         IF (PRESENT(aliases)) aliases = potential%aliases
1038         IF (PRESENT(description)) description = potential%description
1039
1040         IF (PRESENT(elec_conf)) elec_conf => potential%elec_conf
1041
1042         IF (PRESENT(z)) z = potential%z
1043         IF (PRESENT(zeff)) zeff = potential%zeff
1044         IF (PRESENT(zeff_correction)) zeff_correction = potential%zeff_correction
1045         IF (PRESENT(alpha_core_charge)) alpha_core_charge = potential%alpha_core_charge
1046         IF (PRESENT(ccore_charge)) ccore_charge = potential%ccore_charge
1047         IF (PRESENT(core_charge_radius)) core_charge_radius = potential%core_charge_radius
1048
1049         IF (PRESENT(ppl_radius)) ppl_radius = potential%ppl_radius
1050         IF (PRESENT(ppnl_radius)) ppnl_radius = potential%ppnl_radius
1051         IF (PRESENT(ppl_present)) THEN
1052            ppl_present = (potential%nloc > 0 .OR. potential%n_local > 0)
1053         END IF
1054         IF (PRESENT(ppnl_present)) THEN
1055            ppnl_present = ANY(potential%is_nonlocal)
1056         END IF
1057
1058         IF (PRESENT(ecp_local)) ecp_local = potential%ecp_local
1059         IF (PRESENT(n_local)) n_local = potential%n_local
1060         IF (PRESENT(a_local)) a_local => potential%a_local
1061         IF (PRESENT(c_local)) c_local => potential%c_local
1062
1063         IF (PRESENT(nloc)) nloc = potential%nloc
1064         IF (PRESENT(nrloc)) nrloc = potential%nrloc
1065         IF (PRESENT(aloc)) aloc = potential%aloc
1066         IF (PRESENT(bloc)) bloc = potential%bloc
1067
1068         IF (PRESENT(n_nonlocal)) n_nonlocal = potential%n_nonlocal
1069         IF (PRESENT(nppnl)) nppnl = potential%nppnl
1070         IF (PRESENT(lmax)) lmax = potential%lmax
1071         IF (PRESENT(is_nonlocal)) is_nonlocal(:) = potential%is_nonlocal(:)
1072         IF (PRESENT(a_nonlocal)) a_nonlocal => potential%a_nonlocal
1073         IF (PRESENT(c_nonlocal)) c_nonlocal => potential%c_nonlocal
1074         IF (PRESENT(h_nonlocal)) h_nonlocal => potential%h_nonlocal
1075         IF (PRESENT(cprj_ppnl)) cprj_ppnl => potential%cprj_ppnl
1076         IF (PRESENT(vprj_ppnl)) vprj_ppnl => potential%vprj_ppnl
1077
1078         IF (PRESENT(has_nlcc)) has_nlcc = potential%has_nlcc
1079         IF (PRESENT(n_nlcc)) n_nlcc = potential%n_nlcc
1080         IF (PRESENT(a_nlcc)) a_nlcc => potential%a_nlcc
1081         IF (PRESENT(c_nlcc)) c_nlcc => potential%c_nlcc
1082
1083      ELSE
1084
1085         CPABORT("The pointer potential is not associated.")
1086
1087      END IF
1088
1089   END SUBROUTINE get_sgp_potential
1090
1091! **************************************************************************************************
1092!> \brief   Initialise the coefficients of the projectors of the non-local
1093!>          part of the GTH pseudopotential and the transformation matrices
1094!>          for Cartesian overlap integrals between the orbital basis
1095!>          functions and the projector functions.
1096!> \param potential ...
1097!> \date    16.10.2000
1098!> \author  MK
1099!> \version 1.0
1100! **************************************************************************************************
1101   SUBROUTINE init_cprj_ppnl(potential)
1102      TYPE(gth_potential_type), POINTER                  :: potential
1103
1104      CHARACTER(len=*), PARAMETER :: routineN = 'init_cprj_ppnl', routineP = moduleN//':'//routineN
1105
1106      INTEGER                                            :: cpx, cpy, cpz, cx, cy, cz, ico, iprj, &
1107                                                            iprj_ppnl, l, lp, lprj_ppnl, nprj, px, &
1108                                                            py, pz
1109      REAL(KIND=dp)                                      :: alpha_ppnl, cp
1110
1111      nprj = 0
1112
1113      DO l = 0, potential%lppnl
1114         alpha_ppnl = potential%alpha_ppnl(l)
1115         DO iprj_ppnl = 1, potential%nprj_ppnl(l)
1116            lp = iprj_ppnl - 1
1117            lprj_ppnl = l + 2*lp
1118            cp = SQRT(2.0_dp**(2.0_dp*REAL(lprj_ppnl, dp) + 3.5_dp)* &
1119                      alpha_ppnl**(REAL(lprj_ppnl, dp) + 1.5_dp)/ &
1120                      (rootpi*dfac(2*lprj_ppnl + 1)))
1121            potential%cprj_ppnl(iprj_ppnl, l) = cp
1122            DO cx = 0, l
1123               DO cy = 0, l - cx
1124                  cz = l - cx - cy
1125                  iprj = nprj + co(cx, cy, cz)
1126                  DO px = 0, lp
1127                     DO py = 0, lp - px
1128                        pz = lp - px - py
1129                        cpx = cx + 2*px
1130                        cpy = cy + 2*py
1131                        cpz = cz + 2*pz
1132                        ico = coset(cpx, cpy, cpz)
1133                        potential%cprj(ico, iprj) = cp*fac(lp)/(fac(px)*fac(py)*fac(pz))
1134                     END DO
1135                  END DO
1136               END DO
1137            END DO
1138            nprj = nprj + nco(l)
1139         END DO
1140      END DO
1141
1142   END SUBROUTINE init_cprj_ppnl
1143
1144! **************************************************************************************************
1145!> \brief   Initialise a GTH potential data set structure.
1146!> \param potential ...
1147!> \date    27.10.2000
1148!> \author  MK
1149!> \version 1.0
1150! **************************************************************************************************
1151   SUBROUTINE init_gth_potential(potential)
1152      TYPE(gth_potential_type), POINTER                  :: potential
1153
1154      CHARACTER(len=*), PARAMETER :: routineN = 'init_gth_potential', &
1155         routineP = moduleN//':'//routineN
1156
1157      IF (.NOT. ASSOCIATED(potential)) RETURN
1158
1159      IF (potential%nppnl > 0) THEN
1160
1161         !     *** Initialise the projector coefficients of the    ***
1162         !     *** non-local part of the GTH pseudopotential and   ***
1163         !     *** the transformation matrices "pgf" -> "prj_ppnl" ***
1164
1165         CALL init_cprj_ppnl(potential)
1166
1167         !     *** Initialise the h(i,j) projector coefficients of ***
1168         !     *** the non-local part of the GTH pseudopotential   ***
1169
1170         CALL init_vprj_ppnl(potential)
1171
1172      END IF
1173
1174   END SUBROUTINE init_gth_potential
1175
1176! **************************************************************************************************
1177!> \brief   Initialise the h(i,j) projector coefficients of the non-local part
1178!>          of the GTH pseudopotential.
1179!> \param potential ...
1180!> \date    24.10.2000
1181!> \author  MK
1182!> \version 1.0
1183! **************************************************************************************************
1184   SUBROUTINE init_vprj_ppnl(potential)
1185      TYPE(gth_potential_type), POINTER                  :: potential
1186
1187      CHARACTER(len=*), PARAMETER :: routineN = 'init_vprj_ppnl', routineP = moduleN//':'//routineN
1188
1189      INTEGER                                            :: i, ico, iprj, iprj_ppnl, iso, j, jco, &
1190                                                            jprj, jprj_ppnl, l, nprj
1191
1192      nprj = 0
1193
1194      DO l = 0, potential%lppnl
1195         DO iprj_ppnl = 1, potential%nprj_ppnl(l)
1196            iprj = nprj + (iprj_ppnl - 1)*nco(l)
1197            DO jprj_ppnl = 1, potential%nprj_ppnl(l)
1198               jprj = nprj + (jprj_ppnl - 1)*nco(l)
1199               DO ico = 1, nco(l)
1200                  i = iprj + ico
1201                  DO jco = 1, nco(l)
1202                     j = jprj + jco
1203                     DO iso = 1, nso(l)
1204                        potential%vprj_ppnl(i, j) = potential%vprj_ppnl(i, j) + &
1205                                                    orbtramat(l)%slm(iso, ico)* &
1206                                                    potential%hprj_ppnl(iprj_ppnl, &
1207                                                                        jprj_ppnl, l)* &
1208                                                    orbtramat(l)%slm(iso, jco)
1209                     END DO
1210                  END DO
1211               END DO
1212            END DO
1213         END DO
1214         nprj = nprj + potential%nprj_ppnl(l)*nco(l)
1215      END DO
1216
1217   END SUBROUTINE init_vprj_ppnl
1218
1219! **************************************************************************************************
1220!> \brief ...
1221!> \param potential ...
1222!> \param itype ...
1223!> \param zeff ...
1224!> \param zeff_correction ...
1225! **************************************************************************************************
1226   SUBROUTINE init_all_potential(potential, itype, zeff, zeff_correction)
1227
1228      TYPE(all_potential_type), POINTER                  :: potential
1229      CHARACTER(LEN=*), OPTIONAL                         :: itype
1230      REAL(KIND=dp), OPTIONAL                            :: zeff, zeff_correction
1231
1232      CHARACTER(len=*), PARAMETER :: routineN = 'init_all_potential', &
1233         routineP = moduleN//':'//routineN
1234
1235      INTEGER                                            :: dz
1236
1237      IF (.NOT. ASSOCIATED(potential)) RETURN
1238
1239      IF (PRESENT(zeff)) potential%zeff = zeff
1240      IF (PRESENT(zeff_correction)) potential%zeff_correction = zeff_correction
1241      dz = potential%z - INT(potential%zeff - potential%zeff_correction)
1242      SELECT CASE (dz)
1243      CASE DEFAULT
1244      CASE (2)
1245         potential%elec_conf(0) = potential%elec_conf(0) - 2
1246      CASE (10)
1247         potential%elec_conf(0) = potential%elec_conf(0) - 4
1248         potential%elec_conf(1) = potential%elec_conf(1) - 6
1249      CASE (18)
1250         potential%elec_conf(0) = potential%elec_conf(0) - 6
1251         potential%elec_conf(1) = potential%elec_conf(1) - 12
1252      CASE (28)
1253         potential%elec_conf(0) = potential%elec_conf(0) - 6
1254         potential%elec_conf(1) = potential%elec_conf(1) - 12
1255         potential%elec_conf(2) = potential%elec_conf(2) - 10
1256      CASE (30)
1257         potential%elec_conf(0) = potential%elec_conf(0) - 8
1258         potential%elec_conf(1) = potential%elec_conf(1) - 12
1259         potential%elec_conf(2) = potential%elec_conf(2) - 10
1260      CASE (36)
1261         potential%elec_conf(0) = potential%elec_conf(0) - 8
1262         potential%elec_conf(1) = potential%elec_conf(1) - 18
1263         potential%elec_conf(2) = potential%elec_conf(2) - 10
1264      CASE (46)
1265         potential%elec_conf(0) = potential%elec_conf(0) - 8
1266         potential%elec_conf(1) = potential%elec_conf(1) - 18
1267         potential%elec_conf(2) = potential%elec_conf(2) - 20
1268      CASE (48)
1269         potential%elec_conf(0) = potential%elec_conf(0) - 10
1270         potential%elec_conf(1) = potential%elec_conf(1) - 18
1271         potential%elec_conf(2) = potential%elec_conf(2) - 20
1272      CASE (54)
1273         potential%elec_conf(0) = potential%elec_conf(0) - 10
1274         potential%elec_conf(1) = potential%elec_conf(1) - 24
1275         potential%elec_conf(2) = potential%elec_conf(2) - 20
1276      CASE (68)
1277         potential%elec_conf(0) = potential%elec_conf(0) - 10
1278         potential%elec_conf(1) = potential%elec_conf(1) - 24
1279         potential%elec_conf(2) = potential%elec_conf(2) - 20
1280         potential%elec_conf(3) = potential%elec_conf(3) - 14
1281      CASE (78)
1282         potential%elec_conf(0) = potential%elec_conf(0) - 10
1283         potential%elec_conf(1) = potential%elec_conf(1) - 24
1284         potential%elec_conf(2) = potential%elec_conf(2) - 30
1285         potential%elec_conf(3) = potential%elec_conf(3) - 14
1286      CASE (80)
1287         potential%elec_conf(0) = potential%elec_conf(0) - 12
1288         potential%elec_conf(1) = potential%elec_conf(1) - 24
1289         potential%elec_conf(2) = potential%elec_conf(2) - 30
1290         potential%elec_conf(3) = potential%elec_conf(3) - 14
1291      CASE (86)
1292         potential%elec_conf(0) = potential%elec_conf(0) - 12
1293         potential%elec_conf(1) = potential%elec_conf(1) - 30
1294         potential%elec_conf(2) = potential%elec_conf(2) - 30
1295         potential%elec_conf(3) = potential%elec_conf(3) - 14
1296      CASE (100)
1297         potential%elec_conf(0) = potential%elec_conf(0) - 12
1298         potential%elec_conf(1) = potential%elec_conf(1) - 30
1299         potential%elec_conf(2) = potential%elec_conf(2) - 30
1300         potential%elec_conf(3) = potential%elec_conf(3) - 28
1301      END SELECT
1302
1303      IF (PRESENT(itype)) THEN
1304         IF (itype == "BARE") THEN
1305            potential%description(1) = "Bare Coulomb Potential"
1306            IF (dz > 0) THEN
1307               potential%description(2) = "Valence charge only"
1308            ELSE
1309               potential%description(2) = "Full atomic charge"
1310            END IF
1311         ENDIF
1312      END IF
1313
1314   END SUBROUTINE init_all_potential
1315! **************************************************************************************************
1316!> \brief   Initialise a SGP potential data set structure.
1317!> \param potential ...
1318!> \version 1.0
1319! **************************************************************************************************
1320   SUBROUTINE init_sgp_potential(potential)
1321      TYPE(sgp_potential_type), POINTER                  :: potential
1322
1323      CHARACTER(len=*), PARAMETER :: routineN = 'init_sgp_potential', &
1324         routineP = moduleN//':'//routineN
1325
1326      INTEGER                                            :: i1, i2, j1, j2, l, la, lb, n1, n2, nnl, &
1327                                                            nprj
1328      INTEGER, ALLOCATABLE, DIMENSION(:, :)              :: ind1, ind2
1329      REAL(KIND=dp), DIMENSION(:, :), POINTER            :: cprj, hnl
1330      REAL(KIND=dp), DIMENSION(:, :, :), POINTER         :: cn
1331
1332      IF (ASSOCIATED(potential)) THEN
1333         IF (potential%nppnl > 0) THEN
1334            !
1335            IF (ASSOCIATED(potential%cprj_ppnl)) THEN
1336               DEALLOCATE (potential%cprj_ppnl)
1337            END IF
1338            nprj = 0
1339            DO l = 0, potential%lmax
1340               nnl = potential%n_nonlocal
1341               IF (potential%is_nonlocal(l)) nprj = nprj + nnl*nso(l)
1342            END DO
1343            ALLOCATE (potential%cprj_ppnl(potential%nppnl, nprj))
1344            cprj => potential%cprj_ppnl
1345            cprj = 0.0_dp
1346            cn => potential%c_nonlocal
1347            nnl = potential%n_nonlocal
1348            !
1349            n1 = 0
1350            DO la = 0, potential%lmax
1351               n1 = n1 + nnl*nco(la)
1352            END DO
1353            ALLOCATE (ind1(n1, 3))
1354            n1 = 0
1355            DO i1 = 1, nnl
1356               DO la = 0, potential%lmax
1357                  DO j1 = 1, nco(la)
1358                     n1 = n1 + 1
1359                     ind1(n1, 1) = la
1360                     ind1(n1, 2) = j1
1361                     ind1(n1, 3) = i1
1362                  END DO
1363               END DO
1364            END DO
1365            !
1366            n2 = 0
1367            DO i2 = 1, nnl
1368               DO lb = 0, potential%lmax
1369                  IF (.NOT. potential%is_nonlocal(lb)) CYCLE
1370                  n2 = n2 + nso(lb)
1371               END DO
1372            END DO
1373            ALLOCATE (ind2(n2, 3))
1374            n2 = 0
1375            DO i2 = 1, nnl
1376               DO lb = 0, potential%lmax
1377                  IF (.NOT. potential%is_nonlocal(lb)) CYCLE
1378                  DO j2 = 1, nso(lb)
1379                     n2 = n2 + 1
1380                     ind2(n2, 1) = lb
1381                     ind2(n2, 2) = j2
1382                     ind2(n2, 3) = i2
1383                  END DO
1384               END DO
1385            END DO
1386            !
1387            DO n1 = 1, SIZE(ind1, 1)
1388               la = ind1(n1, 1)
1389               j1 = ind1(n1, 2)
1390               i1 = ind1(n1, 3)
1391               DO n2 = 1, SIZE(ind2, 1)
1392                  lb = ind2(n2, 1)
1393                  IF (la /= lb) CYCLE
1394                  j2 = ind2(n2, 2)
1395                  i2 = ind2(n2, 3)
1396                  cprj(n1, n2) = orbtramat(la)%c2s(j2, j1)*cn(i1, i2, la)
1397               END DO
1398            END DO
1399            !
1400            hnl => potential%h_nonlocal
1401            IF (ASSOCIATED(potential%vprj_ppnl)) THEN
1402               DEALLOCATE (potential%vprj_ppnl)
1403            END IF
1404            ALLOCATE (potential%vprj_ppnl(nprj))
1405            potential%vprj_ppnl = 0.0_dp
1406            DO n2 = 1, SIZE(ind2, 1)
1407               lb = ind2(n2, 1)
1408               i2 = ind2(n2, 3)
1409               potential%vprj_ppnl(n2) = hnl(i2, lb)
1410            END DO
1411            !
1412            DEALLOCATE (ind1, ind2)
1413         END IF
1414      END IF
1415
1416   END SUBROUTINE init_sgp_potential
1417
1418! **************************************************************************************************
1419!> \brief   Read an atomic all-electron potential data set.
1420!> \param element_symbol ...
1421!> \param potential_name ...
1422!> \param potential ...
1423!> \param zeff_correction ...
1424!> \param para_env ...
1425!> \param potential_file_name ...
1426!> \param potential_section ...
1427!> \param update_input ...
1428!> \date    14.05.2000
1429!> \author  MK
1430!> \version 1.0
1431! **************************************************************************************************
1432   SUBROUTINE read_all_potential(element_symbol, potential_name, potential, zeff_correction, &
1433                                 para_env, potential_file_name, potential_section, update_input)
1434
1435      CHARACTER(LEN=*), INTENT(IN)                       :: element_symbol, potential_name
1436      TYPE(all_potential_type), POINTER                  :: potential
1437      REAL(KIND=dp)                                      :: zeff_correction
1438      TYPE(cp_para_env_type), POINTER                    :: para_env
1439      CHARACTER(len=default_path_length), INTENT(IN)     :: potential_file_name
1440      TYPE(section_vals_type), POINTER                   :: potential_section
1441      LOGICAL, INTENT(IN)                                :: update_input
1442
1443      CHARACTER(len=*), PARAMETER :: routineN = 'read_all_potential', &
1444         routineP = moduleN//':'//routineN
1445
1446      CHARACTER(LEN=240)                                 :: line
1447      CHARACTER(LEN=242)                                 :: line2
1448      CHARACTER(len=5*default_string_length)             :: line_att
1449      CHARACTER(LEN=LEN(element_symbol))                 :: symbol
1450      CHARACTER(LEN=LEN(element_symbol)+2)               :: symbol2
1451      CHARACTER(LEN=LEN(potential_name))                 :: apname
1452      CHARACTER(LEN=LEN(potential_name)+2)               :: apname2
1453      INTEGER                                            :: irep, l, strlen1, strlen2
1454      INTEGER, DIMENSION(:), POINTER                     :: elec_conf
1455      LOGICAL                                            :: found, is_ok, match, read_from_input
1456      REAL(KIND=dp)                                      :: alpha, r
1457      TYPE(cp_parser_type), POINTER                      :: parser
1458      TYPE(cp_sll_val_type), POINTER                     :: list
1459      TYPE(val_type), POINTER                            :: val
1460
1461      line2 = ""
1462      symbol2 = ""
1463      apname2 = ""
1464      NULLIFY (parser)
1465      CALL cite_reference(Krack2000)
1466
1467      potential%name = potential_name
1468      read_from_input = .FALSE.
1469      CALL section_vals_get(potential_section, explicit=read_from_input)
1470      IF (.NOT. read_from_input) THEN
1471         CALL parser_create(parser, potential_file_name, para_env=para_env)
1472      END IF
1473
1474      !   *** Search for the requested potential in the potential file   ***
1475      !   *** until the potential is found or the end of file is reached ***
1476
1477      apname = potential_name
1478      symbol = element_symbol
1479      irep = 0
1480      search_loop: DO
1481         IF (read_from_input) THEN
1482            NULLIFY (list, val)
1483            found = .TRUE.
1484            CALL section_vals_list_get(potential_section, "_DEFAULT_KEYWORD_", list=list)
1485         ELSE
1486            CALL parser_search_string(parser, TRIM(apname), .TRUE., found, line)
1487         END IF
1488         IF (found) THEN
1489            CALL uppercase(symbol)
1490            CALL uppercase(apname)
1491
1492            IF (read_from_input) THEN
1493               match = .TRUE.
1494            ELSE
1495               ! Check both the element symbol and the atomic potential name
1496               match = .FALSE.
1497               CALL uppercase(line)
1498               line2 = " "//line//" "
1499               symbol2 = " "//TRIM(symbol)//" "
1500               apname2 = " "//TRIM(apname)//" "
1501               strlen1 = LEN_TRIM(symbol2) + 1
1502               strlen2 = LEN_TRIM(apname2) + 1
1503
1504               IF ((INDEX(line2, symbol2(:strlen1)) > 0) .AND. &
1505                   (INDEX(line2, apname2(:strlen2)) > 0)) match = .TRUE.
1506            END IF
1507            IF (match) THEN
1508               ! Read the electronic configuration
1509               NULLIFY (elec_conf)
1510               l = 0
1511               CALL reallocate(elec_conf, 0, l)
1512               IF (read_from_input) THEN
1513                  is_ok = cp_sll_val_next(list, val)
1514                  IF (.NOT. is_ok) &
1515                     CALL cp_abort(__LOCATION__, &
1516                                   "Error reading the Potential from input file!!")
1517                  CALL val_get(val, c_val=line_att)
1518                  READ (line_att, *) elec_conf(l)
1519                  CALL remove_word(line_att)
1520                  DO WHILE (LEN_TRIM(line_att) /= 0)
1521                     l = l + 1
1522                     CALL reallocate(elec_conf, 0, l)
1523                     READ (line_att, *) elec_conf(l)
1524                     CALL remove_word(line_att)
1525                  END DO
1526               ELSE
1527                  CALL parser_get_object(parser, elec_conf(l), newline=.TRUE.)
1528                  DO WHILE (parser_test_next_token(parser) == "INT")
1529                     l = l + 1
1530                     CALL reallocate(elec_conf, 0, l)
1531                     CALL parser_get_object(parser, elec_conf(l))
1532                  END DO
1533                  irep = irep + 1
1534                  IF (update_input) THEN
1535                     WRITE (line_att, '(100(1X,I0))') elec_conf
1536                     CALL section_vals_val_set(potential_section, "_DEFAULT_KEYWORD_", i_rep_val=irep, &
1537                                               c_val=TRIM(line_att))
1538                  END IF
1539               END IF
1540
1541               CALL reallocate(potential%elec_conf, 0, l)
1542               potential%elec_conf(:) = elec_conf(:)
1543
1544               potential%zeff_correction = zeff_correction
1545               potential%zeff = REAL(SUM(elec_conf), dp) + zeff_correction
1546
1547               DEALLOCATE (elec_conf)
1548
1549               ! Read r(loc) to define the exponent of the core charge
1550               ! distribution and calculate the corresponding coefficient
1551
1552               IF (read_from_input) THEN
1553                  is_ok = cp_sll_val_next(list, val)
1554                  IF (.NOT. is_ok) &
1555                     CALL cp_abort(__LOCATION__, &
1556                                   "Error reading the Potential from input file!!")
1557                  CALL val_get(val, c_val=line_att)
1558                  READ (line_att, *) r
1559               ELSE
1560                  CALL parser_get_object(parser, r, newline=.TRUE.)
1561                  irep = irep + 1
1562                  IF (update_input) THEN
1563                     WRITE (line_att, '(E24.16)') r
1564                     CALL section_vals_val_set(potential_section, "_DEFAULT_KEYWORD_", i_rep_val=irep, &
1565                                               c_val=TRIM(line_att))
1566                  END IF
1567               END IF
1568               alpha = 1.0_dp/(2.0_dp*r**2)
1569
1570               potential%alpha_core_charge = alpha
1571               potential%ccore_charge = potential%zeff*SQRT((alpha/pi)**3)
1572
1573               EXIT search_loop
1574            END IF
1575         ELSE
1576            ! Stop program, if the end of file is reached
1577            CALL cp_abort(__LOCATION__, &
1578                          "The requested atomic potential <"// &
1579                          TRIM(potential_name)// &
1580                          "> for element <"// &
1581                          TRIM(symbol)// &
1582                          "> was not found in the potential file <"// &
1583                          TRIM(potential_file_name)//">")
1584         END IF
1585      END DO search_loop
1586
1587      IF (.NOT. read_from_input) THEN
1588         ! Dump the potential info the in potential section
1589         IF (match .AND. update_input) THEN
1590            irep = irep + 1
1591            WRITE (line_att, '(A)') "         # Potential name: "//apname2(:strlen2)//" for symbol: "//symbol2(:strlen1)
1592            CALL section_vals_val_set(potential_section, "_DEFAULT_KEYWORD_", i_rep_val=irep, &
1593                                      c_val=TRIM(line_att))
1594            irep = irep + 1
1595            WRITE (line_att, '(A)') "         # Potential read from the potential filename: "//TRIM(potential_file_name)
1596            CALL section_vals_val_set(potential_section, "_DEFAULT_KEYWORD_", i_rep_val=irep, &
1597                                      c_val=TRIM(line_att))
1598         END IF
1599         CALL parser_release(parser)
1600      END IF
1601   END SUBROUTINE read_all_potential
1602
1603! **************************************************************************************************
1604!> \brief   Read an atomic local potential data set.
1605!> \param element_symbol ...
1606!> \param potential_name ...
1607!> \param potential ...
1608!> \param para_env ...
1609!> \param potential_file_name ...
1610!> \param potential_section ...
1611!> \param update_input ...
1612!> \date    24.12.2014
1613!> \author  JGH
1614!> \version 1.0
1615! **************************************************************************************************
1616   SUBROUTINE read_local_potential(element_symbol, potential_name, potential, &
1617                                   para_env, potential_file_name, potential_section, update_input)
1618
1619      CHARACTER(LEN=*), INTENT(IN)                       :: element_symbol, potential_name
1620      TYPE(local_potential_type), POINTER                :: potential
1621      TYPE(cp_para_env_type), POINTER                    :: para_env
1622      CHARACTER(len=default_path_length), INTENT(IN)     :: potential_file_name
1623      TYPE(section_vals_type), POINTER                   :: potential_section
1624      LOGICAL, INTENT(IN)                                :: update_input
1625
1626      CHARACTER(len=*), PARAMETER :: routineN = 'read_local_potential', &
1627         routineP = moduleN//':'//routineN
1628      REAL(KIND=dp), PARAMETER                           :: eps_tpot = 1.e-10_dp
1629
1630      CHARACTER(LEN=240)                                 :: line
1631      CHARACTER(LEN=242)                                 :: line2
1632      CHARACTER(len=5*default_string_length)             :: line_att
1633      CHARACTER(LEN=LEN(element_symbol))                 :: symbol
1634      CHARACTER(LEN=LEN(element_symbol)+2)               :: symbol2
1635      CHARACTER(LEN=LEN(potential_name))                 :: apname
1636      CHARACTER(LEN=LEN(potential_name)+2)               :: apname2
1637      INTEGER                                            :: igau, ipol, irep, l, ngau, npol, &
1638                                                            strlen1, strlen2
1639      LOGICAL                                            :: found, is_ok, match, read_from_input
1640      REAL(KIND=dp), DIMENSION(:), POINTER               :: alpha
1641      REAL(KIND=dp), DIMENSION(:, :), POINTER            :: cval
1642      TYPE(cp_parser_type), POINTER                      :: parser
1643      TYPE(cp_sll_val_type), POINTER                     :: list
1644      TYPE(val_type), POINTER                            :: val
1645
1646      line2 = ""
1647      symbol2 = ""
1648      apname2 = ""
1649      NULLIFY (parser, alpha, cval)
1650
1651      potential%name = potential_name
1652      read_from_input = .FALSE.
1653      CALL section_vals_get(potential_section, explicit=read_from_input)
1654      IF (.NOT. read_from_input) THEN
1655         CALL parser_create(parser, potential_file_name, para_env=para_env)
1656      END IF
1657
1658      !   *** Search for the requested potential in the potential file   ***
1659      !   *** until the potential is found or the end of file is reached ***
1660
1661      apname = potential_name
1662      symbol = element_symbol
1663      irep = 0
1664      search_loop: DO
1665         IF (read_from_input) THEN
1666            NULLIFY (list, val)
1667            found = .TRUE.
1668            CALL section_vals_list_get(potential_section, "_DEFAULT_KEYWORD_", list=list)
1669         ELSE
1670            CALL parser_search_string(parser, TRIM(apname), .TRUE., found, line)
1671         END IF
1672         IF (found) THEN
1673            CALL uppercase(symbol)
1674            CALL uppercase(apname)
1675
1676            IF (read_from_input) THEN
1677               match = .TRUE.
1678            ELSE
1679               ! Check both the element symbol and the atomic potential name
1680               match = .FALSE.
1681               CALL uppercase(line)
1682               line2 = " "//line//" "
1683               symbol2 = " "//TRIM(symbol)//" "
1684               apname2 = " "//TRIM(apname)//" "
1685               strlen1 = LEN_TRIM(symbol2) + 1
1686               strlen2 = LEN_TRIM(apname2) + 1
1687
1688               IF ((INDEX(line2, symbol2(:strlen1)) > 0) .AND. &
1689                   (INDEX(line2, apname2(:strlen2)) > 0)) match = .TRUE.
1690            END IF
1691            IF (match) THEN
1692
1693               ! Read ngau and npol
1694               IF (read_from_input) THEN
1695                  is_ok = cp_sll_val_next(list, val)
1696                  IF (.NOT. is_ok) &
1697                     CALL cp_abort(__LOCATION__, &
1698                                   "Error reading the Potential from input file!!")
1699                  CALL val_get(val, c_val=line_att)
1700                  READ (line_att, *) ngau, npol
1701                  CALL remove_word(line_att)
1702               ELSE
1703                  CALL parser_get_object(parser, ngau, newline=.TRUE.)
1704                  CALL parser_get_object(parser, npol)
1705                  irep = irep + 1
1706                  IF (update_input) THEN
1707                     WRITE (line_att, '(2(1X,I0))') ngau, npol
1708                     CALL section_vals_val_set(potential_section, "_DEFAULT_KEYWORD_", i_rep_val=irep, &
1709                                               c_val=TRIM(line_att))
1710                  END IF
1711               END IF
1712
1713               CALL reallocate(alpha, 1, ngau)
1714               CALL reallocate(cval, 1, ngau, 1, npol)
1715               DO igau = 1, ngau
1716                  IF (read_from_input) THEN
1717                     is_ok = cp_sll_val_next(list, val)
1718                     IF (.NOT. is_ok) &
1719                        CALL cp_abort(__LOCATION__, &
1720                                      "Error reading the Potential from input file!!")
1721                     CALL val_get(val, c_val=line_att)
1722                     READ (line_att, *) alpha(igau), (cval(igau, ipol), ipol=1, npol)
1723                  ELSE
1724                     CALL parser_get_object(parser, alpha(igau), newline=.TRUE.)
1725                     DO ipol = 1, npol
1726                        CALL parser_get_object(parser, cval(igau, ipol), newline=.FALSE.)
1727                     END DO
1728                     irep = irep + 1
1729                     IF (update_input) THEN
1730                        WRITE (line_att, '(8(E24.16))') alpha(igau), (cval(igau, ipol), ipol=1, npol)
1731                        CALL section_vals_val_set(potential_section, "_DEFAULT_KEYWORD_", i_rep_val=irep, &
1732                                                  c_val=TRIM(line_att))
1733                     END IF
1734                  END IF
1735               END DO
1736               alpha = 1.0_dp/(2.0_dp*alpha**2)
1737
1738               potential%ngau = ngau
1739               potential%npol = npol
1740
1741               potential%alpha => alpha
1742               potential%cval => cval
1743
1744               potential%radius = 0.0_dp
1745               DO igau = 1, ngau
1746                  DO ipol = 1, npol
1747                     l = 2*(ipol - 1)
1748                     potential%radius = MAX(potential%radius, &
1749                                            exp_radius(l, alpha(igau), eps_tpot, cval(igau, ipol), &
1750                                                       rlow=potential%radius))
1751                  END DO
1752               END DO
1753
1754               EXIT search_loop
1755            END IF
1756         ELSE
1757            ! Stop program, if the end of file is reached
1758            CALL cp_abort(__LOCATION__, &
1759                          "The requested local atomic potential <"// &
1760                          TRIM(potential_name)// &
1761                          "> for element <"// &
1762                          TRIM(symbol)// &
1763                          "> was not found in the potential file <"// &
1764                          TRIM(potential_file_name)//">")
1765         END IF
1766      END DO search_loop
1767
1768      IF (.NOT. read_from_input) THEN
1769         ! Dump the potential info in the potential section
1770         IF (match .AND. update_input) THEN
1771            irep = irep + 1
1772            WRITE (line_att, '(A)') "         # Potential name: "//apname2(:strlen2)//" for symbol: "//symbol2(:strlen1)
1773            CALL section_vals_val_set(potential_section, "_DEFAULT_KEYWORD_", i_rep_val=irep, &
1774                                      c_val=TRIM(line_att))
1775            irep = irep + 1
1776            WRITE (line_att, '(A)') "         # Potential read from the potential filename: "//TRIM(potential_file_name)
1777            CALL section_vals_val_set(potential_section, "_DEFAULT_KEYWORD_", i_rep_val=irep, &
1778                                      c_val=TRIM(line_att))
1779         END IF
1780         CALL parser_release(parser)
1781      END IF
1782
1783   END SUBROUTINE read_local_potential
1784
1785! **************************************************************************************************
1786!> \brief   Read an atomic GTH potential data set.
1787!> \param element_symbol ...
1788!> \param potential_name ...
1789!> \param potential ...
1790!> \param zeff_correction ...
1791!> \param para_env ...
1792!> \param potential_file_name ...
1793!> \param potential_section ...
1794!> \param update_input ...
1795!> \date    14.05.2000
1796!> \par  Literature
1797!>         - S. Goedecker, M. Teter and J. Hutter,
1798!>                Phys. Rev. B 54, 1703 (1996)
1799!>         - C. Hartwigsen, S. Goedecker and J. Hutter,
1800!>                Phys. Rev. B 58, 3641 (1998)
1801!> \author  MK
1802!> \version 1.0
1803! **************************************************************************************************
1804   SUBROUTINE read_gth_potential(element_symbol, potential_name, potential, zeff_correction, &
1805                                 para_env, potential_file_name, potential_section, update_input)
1806
1807      CHARACTER(LEN=*), INTENT(IN)                       :: element_symbol, potential_name
1808      TYPE(gth_potential_type), POINTER                  :: potential
1809      REAL(KIND=dp), INTENT(IN)                          :: zeff_correction
1810      TYPE(cp_para_env_type), POINTER                    :: para_env
1811      CHARACTER(len=default_path_length), INTENT(IN)     :: potential_file_name
1812      TYPE(section_vals_type), POINTER                   :: potential_section
1813      LOGICAL, INTENT(IN)                                :: update_input
1814
1815      CHARACTER(len=*), PARAMETER :: routineN = 'read_gth_potential', &
1816         routineP = moduleN//':'//routineN
1817
1818      CHARACTER(LEN=240)                                 :: line
1819      CHARACTER(LEN=242)                                 :: line2
1820      CHARACTER(len=5*default_string_length)             :: line_att
1821      CHARACTER(LEN=LEN(element_symbol))                 :: symbol
1822      CHARACTER(LEN=LEN(element_symbol)+2)               :: symbol2
1823      CHARACTER(LEN=LEN(potential_name))                 :: apname
1824      CHARACTER(LEN=LEN(potential_name)+2)               :: apname2
1825      INTEGER                                            :: i, ic, ipot, irep, istr, j, l, lppnl, &
1826                                                            lprj_ppnl_max, maxlppl, n, nppnl, &
1827                                                            nprj_ppnl, nprj_ppnl_max, strlen1, &
1828                                                            strlen2
1829      INTEGER, DIMENSION(:), POINTER                     :: elec_conf
1830      LOGICAL                                            :: found, is_ok, match, read_from_input
1831      REAL(KIND=dp)                                      :: alpha, ci, r, rc2
1832      REAL(KIND=dp), DIMENSION(:), POINTER               :: tmp_vals
1833      REAL(KIND=dp), DIMENSION(:, :, :), POINTER         :: hprj_ppnl
1834      TYPE(cp_parser_type), POINTER                      :: parser
1835      TYPE(cp_sll_val_type), POINTER                     :: list
1836      TYPE(val_type), POINTER                            :: val
1837
1838      line2 = ""
1839      symbol2 = ""
1840      apname2 = ""
1841      NULLIFY (parser, tmp_vals)
1842      CALL cite_reference(Goedecker1996)
1843      CALL cite_reference(Hartwigsen1998)
1844      CALL cite_reference(Krack2005)
1845
1846      potential%name = potential_name
1847      potential%aliases = potential_name
1848      read_from_input = .FALSE.
1849      CALL section_vals_get(potential_section, explicit=read_from_input)
1850      IF (.NOT. read_from_input) THEN
1851         CALL parser_create(parser, potential_file_name, para_env=para_env)
1852      END IF
1853
1854      !initialize extended form
1855      potential%lpotextended = .FALSE.
1856      potential%nexp_lpot = 0
1857      potential%lsdpot = .FALSE.
1858      potential%nexp_lsd = 0
1859      potential%nlcc = .FALSE.
1860      potential%nexp_nlcc = 0
1861
1862      !   *** Search for the requested potential in the potential file   ***
1863      !   *** until the potential is found or the end of file is reached ***
1864
1865      apname = potential_name
1866      symbol = element_symbol
1867      irep = 0
1868      search_loop: DO
1869         IF (read_from_input) THEN
1870            NULLIFY (list, val)
1871            found = .TRUE.
1872            CALL section_vals_list_get(potential_section, "_DEFAULT_KEYWORD_", list=list)
1873         ELSE
1874            CALL parser_search_string(parser, TRIM(apname), .TRUE., found, line)
1875         END IF
1876         IF (found) THEN
1877            CALL uppercase(symbol)
1878            CALL uppercase(apname)
1879            IF (read_from_input) THEN
1880               match = .TRUE.
1881            ELSE
1882               !       *** Check both the element symbol and the atomic potential name ***
1883               match = .FALSE.
1884               CALL uppercase(line)
1885               line2 = " "//line//" "
1886               symbol2 = " "//TRIM(symbol)//" "
1887               apname2 = " "//TRIM(apname)//" "
1888               strlen1 = LEN_TRIM(symbol2) + 1
1889               strlen2 = LEN_TRIM(apname2) + 1
1890
1891               i = INDEX(line2, symbol2(:strlen1))
1892               j = INDEX(line2, apname2(:strlen2))
1893               IF (i > 0 .AND. j > 0) THEN
1894                  match = .TRUE.
1895                  i = i + 1 + INDEX(line2(i + 1:), " ")
1896                  potential%aliases = line2(i:) ! copy all names into aliases field
1897               ENDIF
1898            END IF
1899            IF (match) THEN
1900               !         *** Read the electronic configuration ***
1901               NULLIFY (elec_conf)
1902               l = 0
1903               CALL reallocate(elec_conf, 0, l)
1904               IF (read_from_input) THEN
1905                  is_ok = cp_sll_val_next(list, val)
1906                  IF (.NOT. is_ok) &
1907                     CALL cp_abort(__LOCATION__, &
1908                                   "Error reading the Potential from input file!!")
1909                  CALL val_get(val, c_val=line_att)
1910                  READ (line_att, *) elec_conf(l)
1911                  CALL remove_word(line_att)
1912                  DO WHILE (LEN_TRIM(line_att) /= 0)
1913                     l = l + 1
1914                     CALL reallocate(elec_conf, 0, l)
1915                     READ (line_att, *) elec_conf(l)
1916                     CALL remove_word(line_att)
1917                  END DO
1918               ELSE
1919                  CALL parser_get_object(parser, elec_conf(l), newline=.TRUE.)
1920                  DO WHILE (parser_test_next_token(parser) == "INT")
1921                     l = l + 1
1922                     CALL reallocate(elec_conf, 0, l)
1923                     CALL parser_get_object(parser, elec_conf(l))
1924                  END DO
1925                  irep = irep + 1
1926                  IF (update_input) THEN
1927                     WRITE (line_att, '(100(1X,I0))') elec_conf
1928                     CALL section_vals_val_set(potential_section, "_DEFAULT_KEYWORD_", i_rep_val=irep, &
1929                                               c_val=TRIM(line_att))
1930                  END IF
1931               END IF
1932
1933               CALL reallocate(potential%elec_conf, 0, l)
1934               potential%elec_conf(:) = elec_conf(:)
1935
1936               potential%zeff_correction = zeff_correction
1937               potential%zeff = REAL(SUM(elec_conf), dp) + zeff_correction
1938
1939               DEALLOCATE (elec_conf)
1940
1941               !         *** Read r(loc) to define the exponent of the core charge    ***
1942               !         *** distribution and calculate the corresponding coefficient ***
1943
1944               IF (read_from_input) THEN
1945                  is_ok = cp_sll_val_next(list, val)
1946                  IF (.NOT. is_ok) &
1947                     CALL cp_abort(__LOCATION__, &
1948                                   "Error reading the Potential from input file!!")
1949                  CALL val_get(val, c_val=line_att)
1950                  READ (line_att, *) r
1951                  CALL remove_word(line_att)
1952               ELSE
1953                  line_att = ""
1954                  CALL parser_get_object(parser, r, newline=.TRUE.)
1955                  istr = LEN_TRIM(line_att) + 1
1956                  WRITE (line_att(istr:), '(E24.16)') r
1957               END IF
1958               alpha = 1.0_dp/(2.0_dp*r**2)
1959
1960               potential%alpha_core_charge = alpha
1961               potential%ccore_charge = potential%zeff*SQRT((alpha/pi)**3)
1962
1963               potential%alpha_ppl = alpha
1964               potential%cerf_ppl = potential%zeff*SQRT((alpha/pi)**3)
1965
1966               !         *** Read the parameters for the local part ***
1967               !         *** of the GTH pseudopotential (ppl)       ***
1968
1969               IF (read_from_input) THEN
1970                  READ (line_att, *) n
1971                  CALL remove_word(line_att)
1972               ELSE
1973                  CALL parser_get_object(parser, n)
1974                  istr = LEN_TRIM(line_att) + 1
1975                  WRITE (line_att(istr:), '(1X,I0)') n
1976               END IF
1977               potential%nexp_ppl = n
1978               CALL reallocate(potential%cexp_ppl, 1, n)
1979
1980               DO i = 1, n
1981                  IF (read_from_input) THEN
1982                     READ (line_att, *) ci
1983                     CALL remove_word(line_att)
1984                  ELSE
1985                     CALL parser_get_object(parser, ci)
1986                     istr = LEN_TRIM(line_att) + 1
1987                     WRITE (line_att(istr:), '(E24.16)') ci
1988                  END IF
1989                  rc2 = (2.0_dp*potential%alpha_ppl)
1990                  potential%cexp_ppl(i) = rc2**(i - 1)*ci
1991               END DO
1992
1993               IF (.NOT. read_from_input) THEN
1994                  irep = irep + 1
1995                  IF (update_input) THEN
1996                     CALL section_vals_val_set(potential_section, "_DEFAULT_KEYWORD_", i_rep_val=irep, &
1997                                               c_val=TRIM(line_att))
1998                  END IF
1999                  line_att = ""
2000               ELSE
2001                  IF (LEN_TRIM(line_att) /= 0) THEN
2002                     CALL cp_abort(__LOCATION__, &
2003                                   "Error reading the Potential from input file!!")
2004                  END IF
2005               END IF
2006               maxlppl = 2*(n - 1)
2007
2008               IF (maxlppl > -1) CALL init_orbital_pointers(maxlppl)
2009
2010               ! Read extended form of GTH pseudopotential
2011               ! local potential, NLCC, LSD potential
2012               IF (read_from_input) THEN
2013                  DO
2014                     is_ok = cp_sll_val_next(list, val)
2015                     CPASSERT(is_ok)
2016                     CALL val_get(val, c_val=line_att)
2017                     IF (INDEX(line_att, "LPOT") /= 0) THEN
2018                        potential%lpotextended = .TRUE.
2019                        CALL remove_word(line_att)
2020                        READ (line_att, *) potential%nexp_lpot
2021                        n = potential%nexp_lpot
2022                        maxlppl = 2*(n - 1)
2023                        IF (maxlppl > -1) CALL init_orbital_pointers(maxlppl)
2024                        NULLIFY (potential%alpha_lpot, potential%nct_lpot, potential%cval_lpot)
2025                        CALL reallocate(potential%alpha_lpot, 1, n)
2026                        CALL reallocate(potential%nct_lpot, 1, n)
2027                        CALL reallocate(potential%cval_lpot, 1, 4, 1, n)
2028                        DO ipot = 1, potential%nexp_lpot
2029                           is_ok = cp_sll_val_next(list, val)
2030                           CPASSERT(is_ok)
2031                           CALL val_get(val, c_val=line_att)
2032                           READ (line_att, *) r
2033                           potential%alpha_lpot(ipot) = 0.5_dp/(r*r)
2034                           CALL remove_word(line_att)
2035                           READ (line_att, *) potential%nct_lpot(ipot)
2036                           CALL remove_word(line_att)
2037                           DO ic = 1, potential%nct_lpot(ipot)
2038                              READ (line_att, *) ci
2039                              rc2 = (2._dp*potential%alpha_lpot(ipot))**(ic - 1)
2040                              potential%cval_lpot(ic, ipot) = ci*rc2
2041                              CALL remove_word(line_att)
2042                           END DO
2043                        END DO
2044                     ELSEIF (INDEX(line_att, "NLCC") /= 0) THEN
2045                        potential%nlcc = .TRUE.
2046                        CALL remove_word(line_att)
2047                        READ (line_att, *) potential%nexp_nlcc
2048                        n = potential%nexp_nlcc
2049                        NULLIFY (potential%alpha_nlcc, potential%nct_nlcc, potential%cval_nlcc)
2050                        CALL reallocate(potential%alpha_nlcc, 1, n)
2051                        CALL reallocate(potential%nct_nlcc, 1, n)
2052                        CALL reallocate(potential%cval_nlcc, 1, 4, 1, n)
2053                        DO ipot = 1, potential%nexp_nlcc
2054                           is_ok = cp_sll_val_next(list, val)
2055                           CPASSERT(is_ok)
2056                           CALL val_get(val, c_val=line_att)
2057                           READ (line_att, *) potential%alpha_nlcc(ipot)
2058                           CALL remove_word(line_att)
2059                           READ (line_att, *) potential%nct_nlcc(ipot)
2060                           CALL remove_word(line_att)
2061                           DO ic = 1, potential%nct_nlcc(ipot)
2062                              READ (line_att, *) potential%cval_nlcc(ic, ipot)
2063                              !make it compatible with bigdft style
2064                              potential%cval_nlcc(ic, ipot) = potential%cval_nlcc(ic, ipot)/(4.0_dp*pi)
2065                              CALL remove_word(line_att)
2066                           END DO
2067                        END DO
2068                     ELSEIF (INDEX(line_att, "LSD") /= 0) THEN
2069                        potential%lsdpot = .TRUE.
2070                        CALL remove_word(line_att)
2071                        READ (line_att, *) potential%nexp_lsd
2072                        n = potential%nexp_lsd
2073                        NULLIFY (potential%alpha_lsd, potential%nct_lsd, potential%cval_lsd)
2074                        CALL reallocate(potential%alpha_lsd, 1, n)
2075                        CALL reallocate(potential%nct_lsd, 1, n)
2076                        CALL reallocate(potential%cval_lsd, 1, 4, 1, n)
2077                        DO ipot = 1, potential%nexp_lsd
2078                           is_ok = cp_sll_val_next(list, val)
2079                           CPASSERT(is_ok)
2080                           CALL val_get(val, c_val=line_att)
2081                           READ (line_att, *) r
2082                           potential%alpha_lsd(ipot) = 0.5_dp/(r*r)
2083                           CALL remove_word(line_att)
2084                           READ (line_att, *) potential%nct_lsd(ipot)
2085                           CALL remove_word(line_att)
2086                           DO ic = 1, potential%nct_lsd(ipot)
2087                              READ (line_att, *) ci
2088                              rc2 = (2._dp*potential%alpha_lsd(ipot))**(ic - 1)
2089                              potential%cval_lsd(ic, ipot) = ci*rc2
2090                              CALL remove_word(line_att)
2091                           END DO
2092                        END DO
2093                     ELSE
2094                        EXIT
2095                     END IF
2096                  END DO
2097               ELSE
2098                  DO
2099                     CALL parser_get_next_line(parser, 1)
2100                     IF (parser_test_next_token(parser) == "INT") THEN
2101                        EXIT
2102                     ELSEIF (parser_test_next_token(parser) == "STR") THEN
2103                        CALL parser_get_object(parser, line)
2104                        IF (INDEX(LINE, "LPOT") /= 0) THEN
2105                           ! local potential
2106                           potential%lpotextended = .TRUE.
2107                           CALL parser_get_object(parser, potential%nexp_lpot)
2108                           n = potential%nexp_lpot
2109                           NULLIFY (potential%alpha_lpot, potential%nct_lpot, potential%cval_lpot)
2110                           CALL reallocate(potential%alpha_lpot, 1, n)
2111                           CALL reallocate(potential%nct_lpot, 1, n)
2112                           CALL reallocate(potential%cval_lpot, 1, 4, 1, n)
2113                           ! add to input section
2114                           irep = irep + 1
2115                           IF (update_input) THEN
2116                              WRITE (line_att, '(A,1X,I0)') "LPOT", n
2117                              CALL section_vals_val_set(potential_section, "_DEFAULT_KEYWORD_", i_rep_val=irep, &
2118                                                        c_val=TRIM(line_att))
2119                           END IF
2120                           DO ipot = 1, potential%nexp_lpot
2121                              CALL parser_get_object(parser, r, newline=.TRUE.)
2122                              potential%alpha_lpot(ipot) = 0.5_dp/(r*r)
2123                              CALL parser_get_object(parser, potential%nct_lpot(ipot))
2124                              CALL reallocate(tmp_vals, 1, potential%nct_lpot(ipot))
2125                              DO ic = 1, potential%nct_lpot(ipot)
2126                                 CALL parser_get_object(parser, ci)
2127                                 tmp_vals(ic) = ci
2128                                 rc2 = (2._dp*potential%alpha_lpot(ipot))**(ic - 1)
2129                                 potential%cval_lpot(ic, ipot) = ci*rc2
2130                              END DO
2131                              ! add to input section
2132                              irep = irep + 1
2133                              IF (update_input) THEN
2134                                 WRITE (line_att, '(E24.16,1X,I0,100(1X,E24.16))') r, potential%nct_lpot(ipot), &
2135                                    tmp_vals(1:potential%nct_lpot(ipot))
2136                                 CALL section_vals_val_set(potential_section, "_DEFAULT_KEYWORD_", i_rep_val=irep, &
2137                                                           c_val=TRIM(line_att))
2138                              END IF
2139                           END DO
2140                        ELSEIF (INDEX(LINE, "NLCC") /= 0) THEN
2141                           ! NLCC
2142                           potential%nlcc = .TRUE.
2143                           CALL parser_get_object(parser, potential%nexp_nlcc)
2144                           n = potential%nexp_nlcc
2145                           NULLIFY (potential%alpha_nlcc, potential%nct_nlcc, potential%cval_nlcc)
2146                           CALL reallocate(potential%alpha_nlcc, 1, n)
2147                           CALL reallocate(potential%nct_nlcc, 1, n)
2148                           CALL reallocate(potential%cval_nlcc, 1, 4, 1, n)
2149                           ! add to input section
2150                           WRITE (line_att, '(A,1X,I0)') "NLCC", n
2151                           irep = irep + 1
2152                           CALL section_vals_val_set(potential_section, "_DEFAULT_KEYWORD_", i_rep_val=irep, &
2153                                                     c_val=TRIM(line_att))
2154                           DO ipot = 1, potential%nexp_nlcc
2155                              CALL parser_get_object(parser, potential%alpha_nlcc(ipot), newline=.TRUE.)
2156                              CALL parser_get_object(parser, potential%nct_nlcc(ipot))
2157                              CALL reallocate(tmp_vals, 1, potential%nct_nlcc(ipot))
2158                              DO ic = 1, potential%nct_nlcc(ipot)
2159                                 CALL parser_get_object(parser, potential%cval_nlcc(ic, ipot))
2160                                 tmp_vals(ic) = potential%cval_nlcc(ic, ipot)
2161                                 !make it compatible with bigdft style
2162                                 potential%cval_nlcc(ic, ipot) = potential%cval_nlcc(ic, ipot)/(4.0_dp*pi)
2163                              END DO
2164                              ! add to input section
2165                              irep = irep + 1
2166                              IF (update_input) THEN
2167                                 WRITE (line_att, '(E24.16,1X,I0,100(1X,E24.16))') &
2168                                    potential%alpha_nlcc(ipot), potential%nct_nlcc(ipot), &
2169                                    tmp_vals(1:potential%nct_nlcc(ipot))
2170                                 CALL section_vals_val_set(potential_section, "_DEFAULT_KEYWORD_", i_rep_val=irep, &
2171                                                           c_val=TRIM(line_att))
2172                              END IF
2173                           END DO
2174                        ELSEIF (INDEX(LINE, "LSD") /= 0) THEN
2175                           ! LSD potential
2176                           potential%lsdpot = .TRUE.
2177                           CALL parser_get_object(parser, potential%nexp_lsd)
2178                           n = potential%nexp_lsd
2179                           NULLIFY (potential%alpha_lsd, potential%nct_lsd, potential%cval_lsd)
2180                           CALL reallocate(potential%alpha_lsd, 1, n)
2181                           CALL reallocate(potential%nct_lsd, 1, n)
2182                           CALL reallocate(potential%cval_lsd, 1, 4, 1, n)
2183                           ! add to input section
2184                           irep = irep + 1
2185                           IF (update_input) THEN
2186                              WRITE (line_att, '(A,1X,I0)') "LSD", n
2187                              CALL section_vals_val_set(potential_section, "_DEFAULT_KEYWORD_", i_rep_val=irep, &
2188                                                        c_val=TRIM(line_att))
2189                           END IF
2190                           DO ipot = 1, potential%nexp_lsd
2191                              CALL parser_get_object(parser, r, newline=.TRUE.)
2192                              potential%alpha_lsd(ipot) = 0.5_dp/(r*r)
2193                              CALL parser_get_object(parser, potential%nct_lsd(ipot))
2194                              CALL reallocate(tmp_vals, 1, potential%nct_lsd(ipot))
2195                              DO ic = 1, potential%nct_lsd(ipot)
2196                                 CALL parser_get_object(parser, ci)
2197                                 tmp_vals(ic) = ci
2198                                 rc2 = (2._dp*potential%alpha_lsd(ipot))**(ic - 1)
2199                                 potential%cval_lsd(ic, ipot) = ci*rc2
2200                              END DO
2201                              ! add to input section
2202                              irep = irep + 1
2203                              IF (update_input) THEN
2204                                 WRITE (line_att, '(E24.16,1X,I0,100(1X,E24.16))') r, potential%nct_lsd(ipot), &
2205                                    tmp_vals(1:potential%nct_lsd(ipot))
2206                                 CALL section_vals_val_set(potential_section, "_DEFAULT_KEYWORD_", i_rep_val=irep, &
2207                                                           c_val=TRIM(line_att))
2208                              END IF
2209                           END DO
2210                        ELSE
2211                           CALL cp_abort(__LOCATION__, &
2212                                         "Syntax error for <"// &
2213                                         TRIM(element_symbol)// &
2214                                         "> in the atomic potential <"// &
2215                                         TRIM(potential_name)// &
2216                                         "> potential file <"// &
2217                                         TRIM(potential_file_name)//">: "// &
2218                                         "Expected LPOT/NLCC/LSD keyword, got: <"// &
2219                                         TRIM(line)//">")
2220                        END IF
2221                     ELSE
2222                        CALL parser_get_object(parser, line)
2223                        CALL cp_abort(__LOCATION__, &
2224                                      "Syntax error for <"// &
2225                                      TRIM(element_symbol)// &
2226                                      "> in the atomic potential <"// &
2227                                      TRIM(potential_name)// &
2228                                      "> potential file <"// &
2229                                      TRIM(potential_file_name)//">: "// &
2230                                      "Expected LPOT/NLCC/LSD keyword or INTEGER, got: <"// &
2231                                      TRIM(line)//">")
2232                     END IF
2233                  END DO
2234               END IF
2235               !         *** Read the parameters for the non-local  ***
2236               !         *** part of the GTH pseudopotential (ppnl) ***
2237               IF (read_from_input) THEN
2238                  READ (line_att, *) n
2239                  CALL remove_word(line_att)
2240               ELSE
2241                  CALL parser_get_object(parser, n)
2242                  irep = irep + 1
2243                  IF (update_input) THEN
2244                     WRITE (line_att, '(1X,I0)') n
2245                     CALL section_vals_val_set(potential_section, "_DEFAULT_KEYWORD_", i_rep_val=irep, &
2246                                               c_val=TRIM(line_att))
2247                  END IF
2248               END IF
2249               potential%lppnl = n - 1
2250               potential%nppnl = 0
2251
2252               potential%lprj_ppnl_max = n - 1
2253               potential%nprj_ppnl_max = 0
2254
2255               IF (n > 0) THEN
2256
2257                  lppnl = potential%lppnl
2258                  nppnl = potential%nppnl
2259
2260                  CALL init_orbital_pointers(lppnl)
2261
2262                  NULLIFY (hprj_ppnl)
2263
2264                  ! Load the parameter for n non-local projectors
2265
2266                  CALL reallocate(potential%alpha_ppnl, 0, lppnl)
2267                  CALL reallocate(potential%nprj_ppnl, 0, lppnl)
2268
2269                  lprj_ppnl_max = -1
2270                  nprj_ppnl_max = 0
2271
2272                  DO l = 0, lppnl
2273                     IF (read_from_input) THEN
2274                        is_ok = cp_sll_val_next(list, val)
2275                        IF (.NOT. is_ok) &
2276                           CALL cp_abort(__LOCATION__, &
2277                                         "Error reading the Potential from input file!!")
2278                        CALL val_get(val, c_val=line_att)
2279                        READ (line_att, *) r
2280                        CALL remove_word(line_att)
2281                        READ (line_att, *) nprj_ppnl
2282                        CALL remove_word(line_att)
2283                     ELSE
2284                        line_att = ""
2285                        CALL parser_get_object(parser, r, newline=.TRUE.)
2286                        CALL parser_get_object(parser, nprj_ppnl)
2287                        istr = LEN_TRIM(line_att) + 1
2288                        WRITE (line_att(istr:), '(E24.16,1X,I0)') r, nprj_ppnl
2289                     END IF
2290                     IF (r == 0.0_dp .AND. nprj_ppnl /= 0) THEN
2291                        CALL cp_abort(__LOCATION__, &
2292                                      "An error was detected in the atomic potential <"// &
2293                                      TRIM(potential_name)// &
2294                                      "> potential file <"// &
2295                                      TRIM(potential_file_name)//">")
2296                     END IF
2297                     potential%alpha_ppnl(l) = 0.0_dp
2298                     IF (r /= 0.0_dp .AND. n /= 0) potential%alpha_ppnl(l) = 1.0_dp/(2.0_dp*r**2)
2299                     potential%nprj_ppnl(l) = nprj_ppnl
2300                     nppnl = nppnl + nprj_ppnl*nco(l)
2301                     IF (nprj_ppnl > nprj_ppnl_max) THEN
2302                        nprj_ppnl_max = nprj_ppnl
2303                        CALL reallocate(hprj_ppnl, 1, nprj_ppnl_max, &
2304                                        1, nprj_ppnl_max, &
2305                                        0, lppnl)
2306                     END IF
2307                     DO i = 1, nprj_ppnl
2308                        IF (i == 1) THEN
2309                           IF (read_from_input) THEN
2310                              READ (line_att, *) hprj_ppnl(i, i, l)
2311                              CALL remove_word(line_att)
2312                           ELSE
2313                              CALL parser_get_object(parser, hprj_ppnl(i, i, l))
2314                              istr = LEN_TRIM(line_att) + 1
2315                              WRITE (line_att(istr:), '(E24.16)') hprj_ppnl(i, i, l)
2316                           END IF
2317                        ELSE
2318                           IF (read_from_input) THEN
2319                              IF (LEN_TRIM(line_att) /= 0) &
2320                                 CALL cp_abort(__LOCATION__, &
2321                                               "Error reading the Potential from input file!!")
2322                              is_ok = cp_sll_val_next(list, val)
2323                              IF (.NOT. is_ok) &
2324                                 CALL cp_abort(__LOCATION__, &
2325                                               "Error reading the Potential from input file!!")
2326                              CALL val_get(val, c_val=line_att)
2327                              READ (line_att, *) hprj_ppnl(i, i, l)
2328                              CALL remove_word(line_att)
2329                           ELSE
2330                              irep = irep + 1
2331                              CALL section_vals_val_set(potential_section, "_DEFAULT_KEYWORD_", i_rep_val=irep, &
2332                                                        c_val=TRIM(line_att))
2333                              line_att = ""
2334                              CALL parser_get_object(parser, hprj_ppnl(i, i, l), newline=.TRUE.)
2335                              istr = LEN_TRIM(line_att) + 1
2336                              WRITE (line_att(istr:), '(E24.16)') hprj_ppnl(i, i, l)
2337                           END IF
2338                        END IF
2339                        DO j = i + 1, nprj_ppnl
2340                           IF (read_from_input) THEN
2341                              READ (line_att, *) hprj_ppnl(i, j, l)
2342                              CALL remove_word(line_att)
2343                           ELSE
2344                              CALL parser_get_object(parser, hprj_ppnl(i, j, l))
2345                              istr = LEN_TRIM(line_att) + 1
2346                              WRITE (line_att(istr:), '(E24.16)') hprj_ppnl(i, j, l)
2347                           END IF
2348                        END DO
2349                     END DO
2350                     IF (.NOT. read_from_input) THEN
2351                        irep = irep + 1
2352                        IF (update_input) THEN
2353                           CALL section_vals_val_set(potential_section, "_DEFAULT_KEYWORD_", i_rep_val=irep, &
2354                                                     c_val=TRIM(line_att))
2355                        END IF
2356                        line_att = ""
2357                     ELSE
2358                        IF (LEN_TRIM(line_att) /= 0) THEN
2359                           CALL cp_abort(__LOCATION__, &
2360                                         "Error reading the Potential from input file!!")
2361                        END IF
2362                     END IF
2363                     IF (nprj_ppnl > 1) THEN
2364                        CALL symmetrize_matrix(hprj_ppnl(:, :, l), "upper_to_lower")
2365                     END IF
2366                     lprj_ppnl_max = MAX(lprj_ppnl_max, l + 2*(nprj_ppnl - 1))
2367                  END DO
2368
2369                  potential%nppnl = nppnl
2370                  CALL init_orbital_pointers(lprj_ppnl_max)
2371
2372                  potential%lprj_ppnl_max = lprj_ppnl_max
2373                  potential%nprj_ppnl_max = nprj_ppnl_max
2374                  CALL reallocate(potential%hprj_ppnl, 1, nprj_ppnl_max, &
2375                                  1, nprj_ppnl_max, &
2376                                  0, lppnl)
2377                  potential%hprj_ppnl(:, :, :) = hprj_ppnl(:, :, :)
2378
2379                  CALL reallocate(potential%cprj, 1, ncoset(lprj_ppnl_max), 1, nppnl)
2380                  CALL reallocate(potential%cprj_ppnl, 1, nprj_ppnl_max, 0, lppnl)
2381                  CALL reallocate(potential%vprj_ppnl, 1, nppnl, 1, nppnl)
2382
2383                  DEALLOCATE (hprj_ppnl)
2384               END IF
2385               EXIT search_loop
2386            END IF
2387         ELSE
2388            ! Stop program, if the end of file is reached
2389            CALL cp_abort(__LOCATION__, &
2390                          "The requested atomic potential <"// &
2391                          TRIM(potential_name)// &
2392                          "> for element <"// &
2393                          TRIM(symbol)// &
2394                          "> was not found in the potential file <"// &
2395                          TRIM(potential_file_name)//">")
2396         END IF
2397      END DO search_loop
2398
2399      IF (.NOT. read_from_input) THEN
2400         ! Dump the potential info the in potential section
2401         IF (match .AND. update_input) THEN
2402            irep = irep + 1
2403            WRITE (line_att, '(A)') "         # Potential name: "//apname2(:strlen2)//" for symbol: "//symbol2(:strlen1)
2404            CALL section_vals_val_set(potential_section, "_DEFAULT_KEYWORD_", i_rep_val=irep, &
2405                                      c_val=TRIM(line_att))
2406            irep = irep + 1
2407            WRITE (line_att, '(A)') "         # Potential read from the potential filename: "//TRIM(potential_file_name)
2408            CALL section_vals_val_set(potential_section, "_DEFAULT_KEYWORD_", i_rep_val=irep, &
2409                                      c_val=TRIM(line_att))
2410         END IF
2411         CALL parser_release(parser)
2412      END IF
2413
2414      IF (ASSOCIATED(tmp_vals)) DEALLOCATE (tmp_vals)
2415
2416   END SUBROUTINE read_gth_potential
2417
2418! **************************************************************************************************
2419!> \brief ...
2420!> \param potential ...
2421!> \param z ...
2422!> \param zeff_correction ...
2423! **************************************************************************************************
2424   SUBROUTINE set_default_all_potential(potential, z, zeff_correction)
2425
2426      TYPE(all_potential_type), POINTER                  :: potential
2427      INTEGER, INTENT(IN)                                :: z
2428      REAL(KIND=dp)                                      :: zeff_correction
2429
2430      CHARACTER(len=*), PARAMETER :: routineN = 'set_default_all_potential', &
2431         routineP = moduleN//':'//routineN
2432
2433      CHARACTER(LEN=default_string_length)               :: name
2434      INTEGER, DIMENSION(:), POINTER                     :: elec_conf
2435      REAL(KIND=dp)                                      :: alpha, alpha_core_charge, ccore_charge, &
2436                                                            core_charge_radius, r, zeff
2437
2438      ALLOCATE (elec_conf(0:3))
2439      elec_conf(0:3) = ptable(z)%e_conv(0:3)
2440      zeff = REAL(SUM(elec_conf), dp) + zeff_correction
2441      name = ptable(z)%name
2442
2443      r = ptable(z)%covalent_radius*0.5_dp
2444      r = MAX(r, 0.2_dp)
2445      r = MIN(r, 1.0_dp)
2446      alpha = 1.0_dp/(2.0_dp*r**2)
2447
2448      core_charge_radius = r
2449      alpha_core_charge = alpha
2450      ccore_charge = zeff*SQRT((alpha/pi)**3)
2451
2452      CALL set_all_potential(potential, &
2453                             name=name, &
2454                             alpha_core_charge=alpha_core_charge, &
2455                             ccore_charge=ccore_charge, &
2456                             core_charge_radius=core_charge_radius, &
2457                             z=z, &
2458                             zeff=zeff, &
2459                             zeff_correction=zeff_correction, &
2460                             elec_conf=elec_conf)
2461
2462      DEALLOCATE (elec_conf)
2463
2464   END SUBROUTINE set_default_all_potential
2465
2466! **************************************************************************************************
2467!> \brief   Set the attributes of an all-electron potential data set.
2468!> \param potential ...
2469!> \param name ...
2470!> \param alpha_core_charge ...
2471!> \param ccore_charge ...
2472!> \param core_charge_radius ...
2473!> \param z ...
2474!> \param zeff ...
2475!> \param zeff_correction ...
2476!> \param elec_conf ...
2477!> \date    11.01.2002
2478!> \author  MK
2479!> \version 1.0
2480! **************************************************************************************************
2481   SUBROUTINE set_all_potential(potential, name, alpha_core_charge, &
2482                                ccore_charge, core_charge_radius, z, zeff, &
2483                                zeff_correction, elec_conf)
2484
2485      TYPE(all_potential_type), POINTER                  :: potential
2486      CHARACTER(LEN=default_string_length), INTENT(IN), &
2487         OPTIONAL                                        :: name
2488      REAL(KIND=dp), INTENT(IN), OPTIONAL                :: alpha_core_charge, ccore_charge, &
2489                                                            core_charge_radius
2490      INTEGER, INTENT(IN), OPTIONAL                      :: z
2491      REAL(KIND=dp), INTENT(IN), OPTIONAL                :: zeff, zeff_correction
2492      INTEGER, DIMENSION(:), OPTIONAL, POINTER           :: elec_conf
2493
2494      CHARACTER(len=*), PARAMETER :: routineN = 'set_all_potential', &
2495         routineP = moduleN//':'//routineN
2496
2497      IF (ASSOCIATED(potential)) THEN
2498
2499         IF (PRESENT(name)) potential%name = name
2500         IF (PRESENT(alpha_core_charge)) &
2501            potential%alpha_core_charge = alpha_core_charge
2502         IF (PRESENT(ccore_charge)) potential%ccore_charge = ccore_charge
2503         IF (PRESENT(core_charge_radius)) &
2504            potential%core_charge_radius = core_charge_radius
2505         IF (PRESENT(z)) potential%z = z
2506         IF (PRESENT(zeff)) potential%zeff = zeff
2507         IF (PRESENT(zeff_correction)) potential%zeff_correction = zeff_correction
2508         IF (PRESENT(elec_conf)) THEN
2509            IF (.NOT. ASSOCIATED(potential%elec_conf)) THEN
2510               CALL reallocate(potential%elec_conf, 0, SIZE(elec_conf) - 1)
2511            END IF
2512            potential%elec_conf(:) = elec_conf(:)
2513         END IF
2514
2515      ELSE
2516
2517         CPABORT("The pointer potential is not associated")
2518
2519      END IF
2520
2521   END SUBROUTINE set_all_potential
2522
2523! **************************************************************************************************
2524!> \brief   Set the attributes of an atomic local potential data set.
2525!> \param potential ...
2526!> \param name ...
2527!> \param alpha ...
2528!> \param cval ...
2529!> \param radius ...
2530!> \date    24.01.2014
2531!> \author  JGH
2532!> \version 1.0
2533! **************************************************************************************************
2534   SUBROUTINE set_local_potential(potential, name, alpha, cval, radius)
2535
2536      TYPE(local_potential_type), POINTER                :: potential
2537      CHARACTER(LEN=default_string_length), INTENT(IN), &
2538         OPTIONAL                                        :: name
2539      REAL(KIND=dp), DIMENSION(:), OPTIONAL, POINTER     :: alpha
2540      REAL(KIND=dp), DIMENSION(:, :), OPTIONAL, POINTER  :: cval
2541      REAL(KIND=dp), INTENT(IN), OPTIONAL                :: radius
2542
2543      CHARACTER(len=*), PARAMETER :: routineN = 'set_local_potential', &
2544         routineP = moduleN//':'//routineN
2545
2546      IF (ASSOCIATED(potential)) THEN
2547
2548         IF (PRESENT(name)) potential%name = name
2549         IF (PRESENT(alpha)) potential%alpha => alpha
2550         IF (PRESENT(cval)) potential%cval => cval
2551         IF (PRESENT(radius)) potential%radius = radius
2552
2553      ELSE
2554
2555         CPABORT("The pointer potential is not associated")
2556
2557      END IF
2558
2559   END SUBROUTINE set_local_potential
2560
2561! **************************************************************************************************
2562!> \brief   Set the attributes of an effective charge and inducible point
2563!>          dipole potential data set.
2564!> \param potential ...
2565!> \param apol ...
2566!> \param cpol ...
2567!> \param qeff ...
2568!> \param mm_radius ...
2569!> \param qmmm_corr_radius ...
2570!> \param qmmm_radius ...
2571!> \date    05.03.2010
2572!> \author  Toon.Verstraelen@gmail.com
2573! **************************************************************************************************
2574   SUBROUTINE set_fist_potential(potential, apol, cpol, qeff, mm_radius, &
2575                                 qmmm_corr_radius, qmmm_radius)
2576      TYPE(fist_potential_type), POINTER                 :: potential
2577      REAL(KIND=dp), INTENT(IN), OPTIONAL                :: apol, cpol, qeff, mm_radius, &
2578                                                            qmmm_corr_radius, qmmm_radius
2579
2580      CHARACTER(len=*), PARAMETER :: routineN = 'set_fist_potential', &
2581         routineP = moduleN//':'//routineN
2582
2583      IF (ASSOCIATED(potential)) THEN
2584
2585         IF (PRESENT(apol)) potential%apol = apol
2586         IF (PRESENT(cpol)) potential%cpol = cpol
2587         IF (PRESENT(mm_radius)) potential%mm_radius = mm_radius
2588         IF (PRESENT(qeff)) potential%qeff = qeff
2589         IF (PRESENT(qmmm_corr_radius)) potential%qmmm_corr_radius = qmmm_corr_radius
2590         IF (PRESENT(qmmm_radius)) potential%qmmm_radius = qmmm_radius
2591
2592      ELSE
2593
2594         CPABORT("The pointer potential is not associated")
2595
2596      END IF
2597
2598   END SUBROUTINE set_fist_potential
2599
2600! **************************************************************************************************
2601!> \brief   Set the attributes of a GTH potential data set.
2602!> \param potential ...
2603!> \param name ...
2604!> \param alpha_core_charge ...
2605!> \param alpha_ppl ...
2606!> \param ccore_charge ...
2607!> \param cerf_ppl ...
2608!> \param core_charge_radius ...
2609!> \param ppl_radius ...
2610!> \param ppnl_radius ...
2611!> \param lppnl ...
2612!> \param lprj_ppnl_max ...
2613!> \param nexp_ppl ...
2614!> \param nppnl ...
2615!> \param nprj_ppnl_max ...
2616!> \param z ...
2617!> \param zeff ...
2618!> \param zeff_correction ...
2619!> \param alpha_ppnl ...
2620!> \param cexp_ppl ...
2621!> \param elec_conf ...
2622!> \param nprj_ppnl ...
2623!> \param cprj ...
2624!> \param cprj_ppnl ...
2625!> \param vprj_ppnl ...
2626!> \param hprj_ppnl ...
2627!> \date    11.01.2002
2628!> \author  MK
2629!> \version 1.0
2630! **************************************************************************************************
2631   SUBROUTINE set_gth_potential(potential, name, alpha_core_charge, alpha_ppl, &
2632                                ccore_charge, cerf_ppl, core_charge_radius, &
2633                                ppl_radius, ppnl_radius, lppnl, lprj_ppnl_max, &
2634                                nexp_ppl, nppnl, nprj_ppnl_max, z, zeff, zeff_correction, &
2635                                alpha_ppnl, cexp_ppl, elec_conf, nprj_ppnl, cprj, cprj_ppnl, &
2636                                vprj_ppnl, hprj_ppnl)
2637      TYPE(gth_potential_type), POINTER                  :: potential
2638      CHARACTER(LEN=default_string_length), INTENT(IN), &
2639         OPTIONAL                                        :: name
2640      REAL(KIND=dp), INTENT(IN), OPTIONAL                :: alpha_core_charge, alpha_ppl, &
2641                                                            ccore_charge, cerf_ppl, &
2642                                                            core_charge_radius, ppl_radius, &
2643                                                            ppnl_radius
2644      INTEGER, INTENT(IN), OPTIONAL                      :: lppnl, lprj_ppnl_max, nexp_ppl, nppnl, &
2645                                                            nprj_ppnl_max, z
2646      REAL(KIND=dp), INTENT(IN), OPTIONAL                :: zeff, zeff_correction
2647      REAL(KIND=dp), DIMENSION(:), OPTIONAL, POINTER     :: alpha_ppnl, cexp_ppl
2648      INTEGER, DIMENSION(:), OPTIONAL, POINTER           :: elec_conf, nprj_ppnl
2649      REAL(KIND=dp), DIMENSION(:, :), OPTIONAL, POINTER  :: cprj, cprj_ppnl, vprj_ppnl
2650      REAL(KIND=dp), DIMENSION(:, :, :), OPTIONAL, &
2651         POINTER                                         :: hprj_ppnl
2652
2653      CHARACTER(len=*), PARAMETER :: routineN = 'set_gth_potential', &
2654         routineP = moduleN//':'//routineN
2655
2656      IF (ASSOCIATED(potential)) THEN
2657
2658         IF (PRESENT(name)) potential%name = name
2659         IF (PRESENT(alpha_core_charge)) &
2660            potential%alpha_core_charge = alpha_core_charge
2661         IF (PRESENT(alpha_ppl)) potential%alpha_ppl = alpha_ppl
2662         IF (PRESENT(ccore_charge)) potential%ccore_charge = ccore_charge
2663         IF (PRESENT(cerf_ppl)) potential%cerf_ppl = cerf_ppl
2664         IF (PRESENT(core_charge_radius)) &
2665            potential%core_charge_radius = core_charge_radius
2666         IF (PRESENT(ppl_radius)) potential%ppl_radius = ppl_radius
2667         IF (PRESENT(ppnl_radius)) potential%ppnl_radius = ppnl_radius
2668         IF (PRESENT(lppnl)) potential%lppnl = lppnl
2669         IF (PRESENT(lprj_ppnl_max)) potential%lprj_ppnl_max = lprj_ppnl_max
2670         IF (PRESENT(nexp_ppl)) potential%nexp_ppl = nexp_ppl
2671         IF (PRESENT(nppnl)) potential%nppnl = nppnl
2672         IF (PRESENT(nprj_ppnl_max)) potential%nprj_ppnl_max = nprj_ppnl_max
2673         IF (PRESENT(z)) potential%z = z
2674         IF (PRESENT(zeff)) potential%zeff = zeff
2675         IF (PRESENT(zeff_correction)) potential%zeff_correction = zeff_correction
2676         IF (PRESENT(alpha_ppnl)) potential%alpha_ppnl => alpha_ppnl
2677         IF (PRESENT(cexp_ppl)) potential%cexp_ppl => cexp_ppl
2678         IF (PRESENT(elec_conf)) THEN
2679            IF (ASSOCIATED(potential%elec_conf)) THEN
2680               DEALLOCATE (potential%elec_conf)
2681            ENDIF
2682            ALLOCATE (potential%elec_conf(0:SIZE(elec_conf) - 1))
2683            potential%elec_conf(:) = elec_conf(:)
2684         ENDIF
2685         IF (PRESENT(nprj_ppnl)) potential%nprj_ppnl => nprj_ppnl
2686         IF (PRESENT(cprj)) potential%cprj => cprj
2687         IF (PRESENT(cprj_ppnl)) potential%cprj_ppnl => cprj_ppnl
2688         IF (PRESENT(vprj_ppnl)) potential%vprj_ppnl => vprj_ppnl
2689         IF (PRESENT(hprj_ppnl)) potential%hprj_ppnl => hprj_ppnl
2690
2691      ELSE
2692
2693         CPABORT("The pointer potential is not associated")
2694
2695      END IF
2696
2697   END SUBROUTINE set_gth_potential
2698
2699! **************************************************************************************************
2700!> \brief ...
2701!> \param potential ...
2702!> \param name ...
2703!> \param description ...
2704!> \param aliases ...
2705!> \param elec_conf ...
2706!> \param z ...
2707!> \param zeff ...
2708!> \param zeff_correction ...
2709!> \param alpha_core_charge ...
2710!> \param ccore_charge ...
2711!> \param core_charge_radius ...
2712!> \param ppl_radius ...
2713!> \param ppnl_radius ...
2714!> \param ecp_local ...
2715!> \param n_local ...
2716!> \param a_local ...
2717!> \param c_local ...
2718!> \param nloc ...
2719!> \param nrloc ...
2720!> \param aloc ...
2721!> \param bloc ...
2722!> \param n_nonlocal ...
2723!> \param nppnl ...
2724!> \param lmax ...
2725!> \param is_nonlocal ...
2726!> \param a_nonlocal ...
2727!> \param h_nonlocal ...
2728!> \param c_nonlocal ...
2729!> \param has_nlcc ...
2730!> \param n_nlcc ...
2731!> \param a_nlcc ...
2732!> \param c_nlcc ...
2733! **************************************************************************************************
2734   SUBROUTINE set_sgp_potential(potential, name, description, aliases, elec_conf, &
2735                                z, zeff, zeff_correction, alpha_core_charge, &
2736                                ccore_charge, core_charge_radius, &
2737                                ppl_radius, ppnl_radius, &
2738                                ecp_local, n_local, a_local, c_local, &
2739                                nloc, nrloc, aloc, bloc, &
2740                                n_nonlocal, nppnl, lmax, is_nonlocal, a_nonlocal, h_nonlocal, c_nonlocal, &
2741                                has_nlcc, n_nlcc, a_nlcc, c_nlcc)
2742
2743      TYPE(sgp_potential_type), POINTER                  :: potential
2744      CHARACTER(LEN=default_string_length), INTENT(IN), &
2745         OPTIONAL                                        :: name
2746      CHARACTER(LEN=default_string_length), &
2747         DIMENSION(4), INTENT(IN), OPTIONAL              :: description
2748      CHARACTER(LEN=default_string_length), INTENT(IN), &
2749         OPTIONAL                                        :: aliases
2750      INTEGER, DIMENSION(:), OPTIONAL, POINTER           :: elec_conf
2751      INTEGER, INTENT(IN), OPTIONAL                      :: z
2752      REAL(KIND=dp), INTENT(IN), OPTIONAL                :: zeff, zeff_correction, &
2753                                                            alpha_core_charge, ccore_charge, &
2754                                                            core_charge_radius, ppl_radius, &
2755                                                            ppnl_radius
2756      LOGICAL, INTENT(IN), OPTIONAL                      :: ecp_local
2757      INTEGER, INTENT(IN), OPTIONAL                      :: n_local
2758      REAL(KIND=dp), DIMENSION(:), OPTIONAL, POINTER     :: a_local, c_local
2759      INTEGER, INTENT(IN), OPTIONAL                      :: nloc
2760      INTEGER, DIMENSION(1:10), INTENT(IN), OPTIONAL     :: nrloc
2761      REAL(dp), DIMENSION(1:10), INTENT(IN), OPTIONAL    :: aloc, bloc
2762      INTEGER, INTENT(IN), OPTIONAL                      :: n_nonlocal, nppnl, lmax
2763      LOGICAL, DIMENSION(0:5), INTENT(IN), OPTIONAL      :: is_nonlocal
2764      REAL(KIND=dp), DIMENSION(:), OPTIONAL, POINTER     :: a_nonlocal
2765      REAL(KIND=dp), DIMENSION(:, :), OPTIONAL, POINTER  :: h_nonlocal
2766      REAL(KIND=dp), DIMENSION(:, :, :), OPTIONAL, &
2767         POINTER                                         :: c_nonlocal
2768      LOGICAL, INTENT(IN), OPTIONAL                      :: has_nlcc
2769      INTEGER, INTENT(IN), OPTIONAL                      :: n_nlcc
2770      REAL(KIND=dp), DIMENSION(:), OPTIONAL, POINTER     :: a_nlcc, c_nlcc
2771
2772      CHARACTER(len=*), PARAMETER :: routineN = 'set_sgp_potential', &
2773         routineP = moduleN//':'//routineN
2774
2775      IF (ASSOCIATED(potential)) THEN
2776
2777         IF (PRESENT(name)) potential%name = name
2778         IF (PRESENT(aliases)) potential%aliases = aliases
2779         IF (PRESENT(description)) potential%description = description
2780
2781         IF (PRESENT(elec_conf)) THEN
2782            IF (ASSOCIATED(potential%elec_conf)) THEN
2783               DEALLOCATE (potential%elec_conf)
2784            ENDIF
2785            ALLOCATE (potential%elec_conf(0:SIZE(elec_conf) - 1))
2786            potential%elec_conf(:) = elec_conf(:)
2787         ENDIF
2788
2789         IF (PRESENT(z)) potential%z = z
2790         IF (PRESENT(zeff)) potential%zeff = zeff
2791         IF (PRESENT(zeff_correction)) potential%zeff_correction = zeff_correction
2792         IF (PRESENT(alpha_core_charge)) potential%alpha_core_charge = alpha_core_charge
2793         IF (PRESENT(ccore_charge)) potential%ccore_charge = ccore_charge
2794         IF (PRESENT(core_charge_radius)) potential%core_charge_radius = core_charge_radius
2795
2796         IF (PRESENT(ppl_radius)) potential%ppl_radius = ppl_radius
2797         IF (PRESENT(ppnl_radius)) potential%ppnl_radius = ppnl_radius
2798
2799         IF (PRESENT(ecp_local)) potential%ecp_local = ecp_local
2800         IF (PRESENT(n_local)) potential%n_local = n_local
2801         IF (PRESENT(a_local)) potential%a_local => a_local
2802         IF (PRESENT(c_local)) potential%c_local => c_local
2803
2804         IF (PRESENT(nloc)) potential%nloc = nloc
2805         IF (PRESENT(nrloc)) potential%nrloc = nrloc
2806         IF (PRESENT(aloc)) potential%aloc = aloc
2807         IF (PRESENT(bloc)) potential%bloc = bloc
2808
2809         IF (PRESENT(n_nonlocal)) potential%n_nonlocal = n_nonlocal
2810         IF (PRESENT(nppnl)) potential%nppnl = nppnl
2811         IF (PRESENT(lmax)) potential%lmax = lmax
2812         IF (PRESENT(is_nonlocal)) potential%is_nonlocal(:) = is_nonlocal(:)
2813         IF (PRESENT(a_nonlocal)) potential%a_nonlocal => a_nonlocal
2814         IF (PRESENT(c_nonlocal)) potential%c_nonlocal => c_nonlocal
2815         IF (PRESENT(h_nonlocal)) potential%h_nonlocal => h_nonlocal
2816
2817         IF (PRESENT(has_nlcc)) potential%has_nlcc = has_nlcc
2818         IF (PRESENT(n_nlcc)) potential%n_nlcc = n_nlcc
2819         IF (PRESENT(a_nlcc)) potential%a_nlcc => a_nlcc
2820         IF (PRESENT(c_nlcc)) potential%c_nlcc => c_nlcc
2821
2822      ELSE
2823
2824         CPABORT("The pointer potential is not associated.")
2825
2826      END IF
2827
2828   END SUBROUTINE set_sgp_potential
2829
2830! **************************************************************************************************
2831!> \brief ...
2832!> \param potential ...
2833!> \param output_unit ...
2834! **************************************************************************************************
2835   SUBROUTINE write_all_potential(potential, output_unit)
2836
2837      ! Write an atomic all-electron potential data set to the output unit.
2838
2839      ! - Creation (09.02.2002,MK)
2840
2841      TYPE(all_potential_type), POINTER                  :: potential
2842      INTEGER, INTENT(in)                                :: output_unit
2843
2844      CHARACTER(LEN=20)                                  :: string
2845
2846      IF (output_unit > 0 .AND. ASSOCIATED(potential)) THEN
2847         WRITE (UNIT=output_unit, FMT="(/,T6,A,T41,A40,/)") &
2848            "AE Potential information for", ADJUSTR(TRIM(potential%name))
2849         WRITE (UNIT=output_unit, FMT="(T8,A,T41,A40)") &
2850            "Description: ", TRIM(potential%description(1)), &
2851            "             ", TRIM(potential%description(2))
2852         WRITE (UNIT=output_unit, FMT="(/,T8,A,T69,F12.6)") &
2853            "Gaussian exponent of the core charge distribution: ", &
2854            potential%alpha_core_charge
2855         WRITE (UNIT=string, FMT="(5I4)") potential%elec_conf
2856         WRITE (UNIT=output_unit, FMT="(T8,A,T61,A20)") &
2857            "Electronic configuration (s p d ...):", &
2858            ADJUSTR(TRIM(string))
2859      END IF
2860
2861   END SUBROUTINE write_all_potential
2862
2863! **************************************************************************************************
2864!> \brief ...
2865!> \param potential ...
2866!> \param output_unit ...
2867! **************************************************************************************************
2868   SUBROUTINE write_local_potential(potential, output_unit)
2869
2870      ! Write an atomic local potential data set to the output unit.
2871
2872      ! - Creation (24.01.2014,JGH)
2873
2874      TYPE(local_potential_type), POINTER                :: potential
2875      INTEGER, INTENT(in)                                :: output_unit
2876
2877      INTEGER                                            :: igau, ipol
2878
2879      IF (output_unit > 0 .AND. ASSOCIATED(potential)) THEN
2880         WRITE (UNIT=output_unit, FMT="(/,T6,A,T41,A40)") &
2881            "Local Potential information for", ADJUSTR(TRIM(potential%name))
2882         WRITE (UNIT=output_unit, FMT="(T8,A,T41,A40)") &
2883            "Description: ", TRIM(potential%description(1))
2884         DO igau = 1, potential%ngau
2885            WRITE (UNIT=output_unit, FMT="(T8,A,F12.6,T50,A,4(T68,I2,F10.4))") &
2886               "Exponent: ", potential%alpha(igau), &
2887               "Coefficients: ", (2*ipol - 2, potential%cval(igau, ipol), ipol=1, potential%npol)
2888         END DO
2889      END IF
2890
2891   END SUBROUTINE write_local_potential
2892
2893! **************************************************************************************************
2894!> \brief ...
2895!> \param potential ...
2896!> \param output_unit ...
2897! **************************************************************************************************
2898   SUBROUTINE write_gth_potential(potential, output_unit)
2899
2900      ! Write an atomic GTH potential data set to the output unit.
2901      ! - Creation (09.02.2002,MK)
2902
2903      TYPE(gth_potential_type), POINTER                  :: potential
2904      INTEGER, INTENT(in)                                :: output_unit
2905
2906      CHARACTER(LEN=20)                                  :: string
2907      INTEGER                                            :: i, j, l
2908      REAL(KIND=dp)                                      :: r
2909
2910      IF (output_unit > 0 .AND. ASSOCIATED(potential)) THEN
2911         WRITE (UNIT=output_unit, FMT="(/,T6,A,T41,A40,/)") &
2912            "GTH Potential information for", ADJUSTR(TRIM(potential%name))
2913         WRITE (UNIT=output_unit, FMT="(T8,A,T41,A40)") &
2914            "Description: ", ADJUSTR(TRIM(potential%description(1))), &
2915            "             ", ADJUSTR(TRIM(potential%description(2))), &
2916            "             ", ADJUSTR(TRIM(potential%description(3))), &
2917            "             ", ADJUSTR(TRIM(potential%description(4)))
2918         WRITE (UNIT=output_unit, FMT="(/,T8,A,T69,F12.6)") &
2919            "Gaussian exponent of the core charge distribution: ", &
2920            potential%alpha_core_charge
2921         WRITE (UNIT=string, FMT="(5I4)") potential%elec_conf
2922         WRITE (UNIT=output_unit, FMT="(T8,A,T61,A20)") &
2923            "Electronic configuration (s p d ...):", &
2924            ADJUSTR(TRIM(string))
2925
2926         r = 1.0_dp/SQRT(2.0_dp*potential%alpha_ppl)
2927
2928         WRITE (UNIT=output_unit, FMT="(/,T8,A,/,/,T27,A,/,T21,5F12.6)") &
2929            "Parameters of the local part of the GTH pseudopotential:", &
2930            "rloc        C1          C2          C3          C4", &
2931            r, (potential%cexp_ppl(i)*r**(2*(i - 1)), i=1, potential%nexp_ppl)
2932
2933         IF (potential%lppnl > -1) THEN
2934            WRITE (UNIT=output_unit, FMT="(/,T8,A,/,/,T20,A,/)") &
2935               "Parameters of the non-local part of the GTH pseudopotential:", &
2936               "l      r(l)      h(i,j,l)"
2937            DO l = 0, potential%lppnl
2938               r = SQRT(0.5_dp/potential%alpha_ppnl(l))
2939               WRITE (UNIT=output_unit, FMT="(T19,I2,5F12.6)") &
2940                  l, r, (potential%hprj_ppnl(1, j, l), j=1, potential%nprj_ppnl(l))
2941               DO i = 2, potential%nprj_ppnl(l)
2942                  WRITE (UNIT=output_unit, FMT="(T33,4F12.6)") &
2943                     (potential%hprj_ppnl(i, j, l), j=1, potential%nprj_ppnl(l))
2944               END DO
2945            END DO
2946         END IF
2947      END IF
2948
2949   END SUBROUTINE write_gth_potential
2950
2951! **************************************************************************************************
2952!> \brief ...
2953!> \param potential ...
2954!> \param output_unit ...
2955! **************************************************************************************************
2956   SUBROUTINE write_sgp_potential(potential, output_unit)
2957
2958      TYPE(sgp_potential_type), POINTER                  :: potential
2959      INTEGER, INTENT(in)                                :: output_unit
2960
2961      CHARACTER(LEN=40)                                  :: string
2962      INTEGER                                            :: i, l
2963
2964      IF (output_unit > 0 .AND. ASSOCIATED(potential)) THEN
2965         WRITE (UNIT=output_unit, FMT="(/,T6,A,T41,A40,/)") &
2966            "SGP Potential information for", ADJUSTR(TRIM(potential%name))
2967         WRITE (UNIT=output_unit, FMT="(T8,A,T25,A56)") &
2968            "Description: ", ADJUSTR(TRIM(potential%description(1))), &
2969            "             ", ADJUSTR(TRIM(potential%description(2))), &
2970            "             ", ADJUSTR(TRIM(potential%description(3))), &
2971            "             ", ADJUSTR(TRIM(potential%description(4)))
2972         WRITE (UNIT=output_unit, FMT="(/,T8,A,T69,F12.6)") &
2973            "Gaussian exponent of the core charge distribution: ", &
2974            potential%alpha_core_charge
2975         WRITE (UNIT=string, FMT="(10I4)") potential%elec_conf
2976         WRITE (UNIT=output_unit, FMT="(T8,A,T61,A20)") &
2977            "Electronic configuration (s p d ...):", &
2978            ADJUSTR(TRIM(string))
2979         IF (potential%ecp_local) THEN
2980            IF (potential%nloc > 0) THEN
2981               WRITE (UNIT=output_unit, FMT="(/,T8,'Local pseudopotential')")
2982               WRITE (UNIT=output_unit, FMT="(T20,'r**(n-2)',T50,'Coefficient',T73,'Exponent')")
2983               DO i = 1, potential%nloc
2984                  WRITE (UNIT=output_unit, FMT="(T20,I5,T47,F14.8,T69,F12.6)") &
2985                     potential%nrloc(i), potential%aloc(i), potential%bloc(i)
2986               END DO
2987            END IF
2988         ELSE
2989            IF (potential%n_local > 0) THEN
2990               WRITE (UNIT=output_unit, FMT="(/,T8,'Local pseudopotential')")
2991               WRITE (UNIT=output_unit, FMT="(T8,A,10(T21,6F10.4,/))") &
2992                  'Exponents:', potential%a_local(1:potential%n_local)
2993               WRITE (UNIT=output_unit, FMT="(T8,A,10(T21,6F10.4,/))") &
2994                  'Coefficients:', potential%c_local(1:potential%n_local)
2995            END IF
2996         END IF
2997         ! nonlocal PP
2998         IF (potential%n_nonlocal > 0) THEN
2999            WRITE (UNIT=output_unit, FMT="(/,T8,'Nonlocal pseudopotential')")
3000            WRITE (UNIT=output_unit, FMT="(T8,A,T71,I10)") 'Total number of projectors:', potential%nppnl
3001            WRITE (UNIT=output_unit, FMT="(T8,A,10(T21,6F10.4,/))") &
3002               'Exponents:', potential%a_nonlocal(1:potential%n_nonlocal)
3003            DO l = 0, potential%lmax
3004               WRITE (UNIT=output_unit, FMT="(T8,'Coupling for l=',I4)") l
3005               WRITE (UNIT=output_unit, FMT="(10(T21,6F10.4,/))") &
3006                  potential%h_nonlocal(1:potential%n_nonlocal, l)
3007            END DO
3008         END IF
3009         !
3010         IF (potential%has_nlcc) THEN
3011            WRITE (UNIT=output_unit, FMT="(/,T8,'Nonlinear Core Correction')")
3012            WRITE (UNIT=output_unit, FMT="(T8,A,10(T21,6F10.4,/))") &
3013               'Exponents:', potential%a_nlcc(1:potential%n_nlcc)
3014            WRITE (UNIT=output_unit, FMT="(T8,A,10(T21,6F10.4,/))") &
3015               'Coefficients:', potential%c_nlcc(1:potential%n_nlcc)
3016         END IF
3017      END IF
3018
3019   END SUBROUTINE write_sgp_potential
3020
3021! **************************************************************************************************
3022!> \brief Copy an all_potential_type to a new, unallocated variable
3023!> \param pot_in the input potential to copy
3024!> \param pot_out the newly copied and allocated potential
3025!> \par History: created 12.2019 (A. Bussy)
3026! **************************************************************************************************
3027   SUBROUTINE copy_all_potential(pot_in, pot_out)
3028
3029      TYPE(all_potential_type), INTENT(IN), POINTER      :: pot_in
3030      TYPE(all_potential_type), INTENT(INOUT), POINTER   :: pot_out
3031
3032      CHARACTER(len=*), PARAMETER :: routineN = 'copy_all_potential', &
3033         routineP = moduleN//':'//routineN
3034
3035      CPASSERT(ASSOCIATED(pot_in))
3036      CALL allocate_all_potential(pot_out)
3037
3038      pot_out%name = pot_in%name
3039      pot_out%alpha_core_charge = pot_in%alpha_core_charge
3040      pot_out%ccore_charge = pot_in%ccore_charge
3041      pot_out%core_charge_radius = pot_in%core_charge_radius
3042      pot_out%zeff = pot_in%zeff
3043      pot_out%zeff_correction = pot_in%zeff_correction
3044      pot_out%z = pot_in%z
3045
3046      IF (ASSOCIATED(pot_in%elec_conf)) THEN
3047         ALLOCATE (pot_out%elec_conf(LBOUND(pot_in%elec_conf, 1):UBOUND(pot_in%elec_conf, 1)))
3048         pot_out%elec_conf(:) = pot_in%elec_conf(:)
3049      END IF
3050
3051   END SUBROUTINE copy_all_potential
3052
3053! **************************************************************************************************
3054!> \brief Copy a gth_potential_type to a new, unallocated variable
3055!> \param pot_in the input potential to copy
3056!> \param pot_out the newly copied and allocated potential
3057!> \par History: created 12.2019 (A. Bussy)
3058! **************************************************************************************************
3059   SUBROUTINE copy_gth_potential(pot_in, pot_out)
3060
3061      TYPE(gth_potential_type), INTENT(IN), POINTER      :: pot_in
3062      TYPE(gth_potential_type), INTENT(INOUT), POINTER   :: pot_out
3063
3064      CHARACTER(len=*), PARAMETER :: routineN = 'copy_gth_potential', &
3065         routineP = moduleN//':'//routineN
3066
3067      CPASSERT(ASSOCIATED(pot_in))
3068      CALL allocate_gth_potential(pot_out)
3069
3070      pot_out%name = pot_in%name
3071      pot_out%aliases = pot_in%aliases
3072      pot_out%alpha_core_charge = pot_in%alpha_core_charge
3073      pot_out%alpha_ppl = pot_in%alpha_ppl
3074      pot_out%ccore_charge = pot_in%ccore_charge
3075      pot_out%cerf_ppl = pot_in%cerf_ppl
3076      pot_out%zeff = pot_in%zeff
3077      pot_out%core_charge_radius = pot_in%core_charge_radius
3078      pot_out%ppl_radius = pot_in%ppl_radius
3079      pot_out%ppnl_radius = pot_in%ppnl_radius
3080      pot_out%zeff_correction = pot_in%zeff_correction
3081      pot_out%lppnl = pot_in%lppnl
3082      pot_out%lprj_ppnl_max = pot_in%lprj_ppnl_max
3083      pot_out%nexp_ppl = pot_in%nexp_ppl
3084      pot_out%nppnl = pot_in%nppnl
3085      pot_out%nprj_ppnl_max = pot_in%nprj_ppnl_max
3086      pot_out%z = pot_in%z
3087      pot_out%nlcc = pot_in%nlcc
3088      pot_out%nexp_nlcc = pot_in%nexp_nlcc
3089      pot_out%lsdpot = pot_in%lsdpot
3090      pot_out%nexp_lsd = pot_in%nexp_lsd
3091      pot_out%lpotextended = pot_in%lpotextended
3092      pot_out%nexp_lpot = pot_in%nexp_lpot
3093
3094      IF (ASSOCIATED(pot_in%alpha_ppnl)) THEN
3095         ALLOCATE (pot_out%alpha_ppnl(LBOUND(pot_in%alpha_ppnl, 1):UBOUND(pot_in%alpha_ppnl, 1)))
3096         pot_out%alpha_ppnl(:) = pot_in%alpha_ppnl(:)
3097      END IF
3098      IF (ASSOCIATED(pot_in%cexp_ppl)) THEN
3099         ALLOCATE (pot_out%cexp_ppl(LBOUND(pot_in%cexp_ppl, 1):UBOUND(pot_in%cexp_ppl, 1)))
3100         pot_out%cexp_ppl(:) = pot_in%cexp_ppl(:)
3101      END IF
3102      IF (ASSOCIATED(pot_in%elec_conf)) THEN
3103         ALLOCATE (pot_out%elec_conf(LBOUND(pot_in%elec_conf, 1):UBOUND(pot_in%elec_conf, 1)))
3104         pot_out%elec_conf(:) = pot_in%elec_conf(:)
3105      END IF
3106      IF (ASSOCIATED(pot_in%nprj_ppnl)) THEN
3107         ALLOCATE (pot_out%nprj_ppnl(LBOUND(pot_in%nprj_ppnl, 1):UBOUND(pot_in%nprj_ppnl, 1)))
3108         pot_out%nprj_ppnl(:) = pot_in%nprj_ppnl(:)
3109      END IF
3110      IF (ASSOCIATED(pot_in%cprj)) THEN
3111         ALLOCATE (pot_out%cprj(LBOUND(pot_in%cprj, 1):UBOUND(pot_in%cprj, 1), &
3112                                LBOUND(pot_in%cprj, 2):UBOUND(pot_in%cprj, 2)))
3113         pot_out%cprj(:, :) = pot_in%cprj(:, :)
3114      END IF
3115      IF (ASSOCIATED(pot_in%cprj_ppnl)) THEN
3116         ALLOCATE (pot_out%cprj_ppnl(LBOUND(pot_in%cprj_ppnl, 1):UBOUND(pot_in%cprj_ppnl, 1), &
3117                                     LBOUND(pot_in%cprj_ppnl, 2):UBOUND(pot_in%cprj_ppnl, 2)))
3118         pot_out%cprj_ppnl(:, :) = pot_in%cprj_ppnl(:, :)
3119      END IF
3120      IF (ASSOCIATED(pot_in%vprj_ppnl)) THEN
3121         ALLOCATE (pot_out%vprj_ppnl(LBOUND(pot_in%vprj_ppnl, 1):UBOUND(pot_in%vprj_ppnl, 1), &
3122                                     LBOUND(pot_in%vprj_ppnl, 2):UBOUND(pot_in%vprj_ppnl, 2)))
3123         pot_out%vprj_ppnl(:, :) = pot_in%vprj_ppnl(:, :)
3124      END IF
3125      IF (ASSOCIATED(pot_in%hprj_ppnl)) THEN
3126         ALLOCATE (pot_out%hprj_ppnl(LBOUND(pot_in%hprj_ppnl, 1):UBOUND(pot_in%hprj_ppnl, 1), &
3127                                     LBOUND(pot_in%hprj_ppnl, 2):UBOUND(pot_in%hprj_ppnl, 2), &
3128                                     LBOUND(pot_in%hprj_ppnl, 3):UBOUND(pot_in%hprj_ppnl, 3)))
3129         pot_out%hprj_ppnl(:, :, :) = pot_in%hprj_ppnl(:, :, :)
3130      END IF
3131      IF (ASSOCIATED(pot_in%alpha_nlcc)) THEN
3132         ALLOCATE (pot_out%alpha_nlcc(LBOUND(pot_in%alpha_nlcc, 1):UBOUND(pot_in%alpha_nlcc, 1)))
3133         pot_out%alpha_nlcc(:) = pot_in%alpha_nlcc(:)
3134      END IF
3135      IF (ASSOCIATED(pot_in%nct_nlcc)) THEN
3136         ALLOCATE (pot_out%nct_nlcc(LBOUND(pot_in%nct_nlcc, 1):UBOUND(pot_in%nct_nlcc, 1)))
3137         pot_out%nct_nlcc(:) = pot_in%nct_nlcc(:)
3138      END IF
3139      IF (ASSOCIATED(pot_in%cval_nlcc)) THEN
3140         ALLOCATE (pot_out%cval_nlcc(LBOUND(pot_in%cval_nlcc, 1):UBOUND(pot_in%cval_nlcc, 1), &
3141                                     LBOUND(pot_in%cval_nlcc, 2):UBOUND(pot_in%cval_nlcc, 2)))
3142         pot_out%cval_nlcc(:, :) = pot_in%cval_nlcc(:, :)
3143      END IF
3144      IF (ASSOCIATED(pot_in%alpha_lsd)) THEN
3145         ALLOCATE (pot_out%alpha_lsd(LBOUND(pot_in%alpha_lsd, 1):UBOUND(pot_in%alpha_lsd, 1)))
3146         pot_out%alpha_lsd(:) = pot_in%alpha_lsd(:)
3147      END IF
3148      IF (ASSOCIATED(pot_in%nct_lsd)) THEN
3149         ALLOCATE (pot_out%nct_lsd(LBOUND(pot_in%nct_lsd, 1):UBOUND(pot_in%nct_lsd, 1)))
3150         pot_out%nct_lsd(:) = pot_in%nct_lsd(:)
3151      END IF
3152      IF (ASSOCIATED(pot_in%cval_lsd)) THEN
3153         ALLOCATE (pot_out%cval_lsd(LBOUND(pot_in%cval_lsd, 1):UBOUND(pot_in%cval_lsd, 1), &
3154                                    LBOUND(pot_in%cval_lsd, 2):UBOUND(pot_in%cval_lsd, 2)))
3155         pot_out%cval_lsd(:, :) = pot_in%cval_lsd(:, :)
3156      END IF
3157      IF (ASSOCIATED(pot_in%alpha_lpot)) THEN
3158         ALLOCATE (pot_out%alpha_lpot(LBOUND(pot_in%alpha_lpot, 1):UBOUND(pot_in%alpha_lpot, 1)))
3159         pot_out%alpha_lpot(:) = pot_in%alpha_lpot(:)
3160      END IF
3161      IF (ASSOCIATED(pot_in%nct_lpot)) THEN
3162         ALLOCATE (pot_out%nct_lpot(LBOUND(pot_in%nct_lpot, 1):UBOUND(pot_in%nct_lpot, 1)))
3163         pot_out%nct_lpot(:) = pot_in%nct_lpot(:)
3164      END IF
3165      IF (ASSOCIATED(pot_in%cval_lpot)) THEN
3166         ALLOCATE (pot_out%cval_lpot(LBOUND(pot_in%cval_lpot, 1):UBOUND(pot_in%cval_lpot, 1), &
3167                                     LBOUND(pot_in%cval_lpot, 2):UBOUND(pot_in%cval_lpot, 2)))
3168         pot_out%cval_lpot(:, :) = pot_in%cval_lpot(:, :)
3169      END IF
3170
3171   END SUBROUTINE copy_gth_potential
3172
3173! **************************************************************************************************
3174!> \brief Copy a sgp_potential_type to a new, unallocated variable
3175!> \param pot_in the input potential to copy
3176!> \param pot_out the newly copied and allocated potential
3177!> \par History: created 12.2019 (A. Bussy)
3178! **************************************************************************************************
3179   SUBROUTINE copy_sgp_potential(pot_in, pot_out)
3180
3181      TYPE(sgp_potential_type), INTENT(IN), POINTER      :: pot_in
3182      TYPE(sgp_potential_type), INTENT(INOUT), POINTER   :: pot_out
3183
3184      CHARACTER(len=*), PARAMETER :: routineN = 'copy_sgp_potential', &
3185         routineP = moduleN//':'//routineN
3186
3187      CPASSERT(ASSOCIATED(pot_in))
3188      CALL allocate_sgp_potential(pot_out)
3189
3190      pot_out%name = pot_in%name
3191      pot_out%aliases = pot_in%aliases
3192      pot_out%z = pot_in%z
3193      pot_out%zeff = pot_in%zeff
3194      pot_out%zeff_correction = pot_in%zeff_correction
3195      pot_out%alpha_core_charge = pot_in%alpha_core_charge
3196      pot_out%ccore_charge = pot_in%ccore_charge
3197      pot_out%core_charge_radius = pot_in%core_charge_radius
3198      pot_out%ppl_radius = pot_in%ppl_radius
3199      pot_out%ppnl_radius = pot_in%ppnl_radius
3200      pot_out%ecp_local = pot_in%ecp_local
3201      pot_out%n_local = pot_in%n_local
3202      pot_out%nloc = pot_in%nloc
3203      pot_out%nrloc = pot_in%nrloc
3204      pot_out%aloc = pot_in%aloc
3205      pot_out%bloc = pot_in%bloc
3206      pot_out%n_nonlocal = pot_in%n_nonlocal
3207      pot_out%nppnl = pot_in%nppnl
3208      pot_out%lmax = pot_in%lmax
3209      pot_out%is_nonlocal = pot_in%is_nonlocal
3210      pot_out%has_nlcc = pot_in%has_nlcc
3211      pot_out%n_nlcc = pot_in%n_nlcc
3212
3213      IF (ASSOCIATED(pot_in%elec_conf)) THEN
3214         ALLOCATE (pot_out%elec_conf(LBOUND(pot_in%elec_conf, 1):UBOUND(pot_in%elec_conf, 1)))
3215         pot_out%elec_conf(:) = pot_in%elec_conf(:)
3216      END IF
3217      IF (ASSOCIATED(pot_in%a_local)) THEN
3218         ALLOCATE (pot_out%a_local(LBOUND(pot_in%a_local, 1):UBOUND(pot_in%a_local, 1)))
3219         pot_out%a_local(:) = pot_in%a_local(:)
3220      END IF
3221      IF (ASSOCIATED(pot_in%c_local)) THEN
3222         ALLOCATE (pot_out%c_local(LBOUND(pot_in%c_local, 1):UBOUND(pot_in%c_local, 1)))
3223         pot_out%c_local(:) = pot_in%c_local(:)
3224      END IF
3225      IF (ASSOCIATED(pot_in%a_nonlocal)) THEN
3226         ALLOCATE (pot_out%a_nonlocal(LBOUND(pot_in%a_nonlocal, 1):UBOUND(pot_in%a_nonlocal, 1)))
3227         pot_out%a_nonlocal(:) = pot_in%a_nonlocal(:)
3228      END IF
3229      IF (ASSOCIATED(pot_in%h_nonlocal)) THEN
3230         ALLOCATE (pot_out%h_nonlocal(LBOUND(pot_in%h_nonlocal, 1):UBOUND(pot_in%h_nonlocal, 1), &
3231                                      LBOUND(pot_in%h_nonlocal, 2):UBOUND(pot_in%h_nonlocal, 2)))
3232         pot_out%h_nonlocal(:, :) = pot_in%h_nonlocal(:, :)
3233      END IF
3234      IF (ASSOCIATED(pot_in%c_nonlocal)) THEN
3235         ALLOCATE (pot_out%c_nonlocal(LBOUND(pot_in%c_nonlocal, 1):UBOUND(pot_in%c_nonlocal, 1), &
3236                                      LBOUND(pot_in%c_nonlocal, 2):UBOUND(pot_in%c_nonlocal, 2), &
3237                                      LBOUND(pot_in%c_nonlocal, 3):UBOUND(pot_in%c_nonlocal, 3)))
3238         pot_out%c_nonlocal(:, :, :) = pot_in%c_nonlocal(:, :, :)
3239      END IF
3240      IF (ASSOCIATED(pot_in%cprj_ppnl)) THEN
3241         ALLOCATE (pot_out%cprj_ppnl(LBOUND(pot_in%cprj_ppnl, 1):UBOUND(pot_in%cprj_ppnl, 1), &
3242                                     LBOUND(pot_in%cprj_ppnl, 2):UBOUND(pot_in%cprj_ppnl, 2)))
3243         pot_out%cprj_ppnl(:, :) = pot_in%cprj_ppnl(:, :)
3244      END IF
3245      IF (ASSOCIATED(pot_in%vprj_ppnl)) THEN
3246         ALLOCATE (pot_out%vprj_ppnl(LBOUND(pot_in%vprj_ppnl, 1):UBOUND(pot_in%vprj_ppnl, 1)))
3247         pot_out%vprj_ppnl(:) = pot_in%vprj_ppnl(:)
3248      END IF
3249      IF (ASSOCIATED(pot_in%a_nlcc)) THEN
3250         ALLOCATE (pot_out%a_nlcc(LBOUND(pot_in%a_nlcc, 1):UBOUND(pot_in%a_nlcc, 1)))
3251         pot_out%a_nlcc(:) = pot_in%a_nlcc(:)
3252      END IF
3253      IF (ASSOCIATED(pot_in%c_nlcc)) THEN
3254         ALLOCATE (pot_out%c_nlcc(LBOUND(pot_in%c_nlcc, 1):UBOUND(pot_in%c_nlcc, 1)))
3255         pot_out%c_nlcc(:) = pot_in%c_nlcc(:)
3256      END IF
3257
3258   END SUBROUTINE copy_sgp_potential
3259
3260END MODULE external_potential_types
3261