1!--------------------------------------------------------------------------------------------------!
2!   CP2K: A general program to perform molecular dynamics simulations                              !
3!   Copyright (C) 2000 - 2019  CP2K developers group                                               !
4!--------------------------------------------------------------------------------------------------!
5
6! **************************************************************************************************
7!> \brief Definition of the semi empirical parameter types.
8!> \author JGH (14.08.2004)
9! **************************************************************************************************
10MODULE semi_empirical_types
11   USE basis_set_types,                 ONLY: deallocate_sto_basis_set,&
12                                              sto_basis_set_type
13   USE cp_log_handling,                 ONLY: cp_get_default_logger,&
14                                              cp_logger_type,&
15                                              cp_to_string
16   USE cp_output_handling,              ONLY: cp_p_file,&
17                                              cp_print_key_finished_output,&
18                                              cp_print_key_should_output,&
19                                              cp_print_key_unit_nr
20   USE dg_types,                        ONLY: dg_type
21   USE input_constants,                 ONLY: &
22        do_method_am1, do_method_mndo, do_method_mndod, do_method_pdg, do_method_pm3, &
23        do_method_pm6, do_method_pm6fm, do_method_pnnl, do_method_rm1, do_se_IS_kdso_d, &
24        do_se_IS_slater
25   USE input_section_types,             ONLY: section_vals_type
26   USE kinds,                           ONLY: default_string_length,&
27                                              dp
28   USE multipole_types,                 ONLY: do_multipole_charge,&
29                                              do_multipole_dipole,&
30                                              do_multipole_none,&
31                                              do_multipole_quadrupole
32   USE physcon,                         ONLY: angstrom,&
33                                              evolt,&
34                                              kcalmol
35   USE pw_pool_types,                   ONLY: pw_pool_type
36   USE semi_empirical_expns3_types,     ONLY: semi_empirical_expns3_p_type,&
37                                              semi_empirical_expns3_release
38   USE semi_empirical_mpole_types,      ONLY: semi_empirical_mpole_p_release,&
39                                              semi_empirical_mpole_p_type
40   USE taper_types,                     ONLY: taper_create,&
41                                              taper_release,&
42                                              taper_type
43#include "./base/base_uses.f90"
44
45   IMPLICIT NONE
46
47   PRIVATE
48
49! *** Global parameters ***
50   CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'semi_empirical_types'
51
52! **************************************************************************************************
53!> \brief Semi-empirical type
54! **************************************************************************************************
55   TYPE semi_empirical_type
56      INTEGER                                :: typ
57      INTEGER                                :: nr
58      INTEGER                                :: core_size, atm_int_size
59      CHARACTER(LEN=default_string_length)   :: name
60      LOGICAL                                :: defined, dorb, extended_basis_set
61      LOGICAL                                :: p_orbitals_on_h
62      INTEGER                                :: z
63      REAL(KIND=dp)                          :: zeff
64      INTEGER                                :: natorb
65      REAL(KIND=dp), DIMENSION(:), POINTER :: beta
66      REAL(KIND=dp), DIMENSION(:), POINTER :: sto_exponents
67      REAL(KIND=dp), DIMENSION(:), POINTER :: zn
68      TYPE(sto_basis_set_type), POINTER      :: basis
69      INTEGER                                :: ngauss
70      REAL(KIND=dp)                        :: eheat
71      REAL(KIND=dp)                        :: uss, upp, udd, uff
72      REAL(KIND=dp)                        :: alp
73      REAL(KIND=dp)                        :: eisol
74      REAL(KIND=dp)                        :: ass, asp, app, de, acoul
75      REAL(KIND=dp)                        :: gss, gsp, gpp, gp2
76      REAL(KIND=dp)                        :: gsd, gpd, gdd
77      REAL(KIND=dp)                        :: hsp
78      REAL(KIND=dp)                        :: dd, qq, am, ad, aq
79      REAL(KIND=dp), DIMENSION(2)           :: pre, d
80      REAL(KIND=dp), DIMENSION(4)           :: fn1, fn2, fn3
81      REAL(KIND=dp), DIMENSION(4, 4)         :: bfn1, bfn2, bfn3
82      REAL(KIND=dp)                        :: f0dd, f2dd, f4dd, f0sd, f0pd, f2pd, &
83                                              g1pd, g2sd, g3pd
84      REAL(KIND=dp), DIMENSION(9)          :: ko
85      REAL(KIND=dp), DIMENSION(6)          :: cs
86      REAL(KIND=dp), DIMENSION(52)         :: onec2el
87      ! Specific for PM6 & PM6-FM
88      REAL(KIND=dp), DIMENSION(0:115)      :: xab
89      REAL(KIND=dp), DIMENSION(0:115)      :: aab
90      REAL(KIND=dp)                        :: a, b, c, rho
91      ! One center - two electron integrals
92      REAL(KIND=dp), DIMENSION(:, :), &
93         POINTER                           :: w
94      TYPE(semi_empirical_mpole_p_type), &
95         POINTER, DIMENSION(:)             :: w_mpole
96      ! 1/R^3 residual integral part
97      TYPE(semi_empirical_expns3_p_type), &
98         POINTER, DIMENSION(:)             :: expns3_int
99   END TYPE semi_empirical_type
100
101   TYPE semi_empirical_p_type
102      TYPE(semi_empirical_type), POINTER    :: se_param
103   END TYPE semi_empirical_p_type
104
105! **************************************************************************************************
106!> \brief  Rotation Matrix Type
107!> \author 05.2008 Teodoro Laino [tlaino] - University of Zurich
108! **************************************************************************************************
109   TYPE rotmat_type
110      ! Value of Rotation Matrices
111      REAL(KIND=dp), DIMENSION(3, 3)      :: sp
112      REAL(KIND=dp), DIMENSION(5, 5)      :: sd
113      REAL(KIND=dp), DIMENSION(6, 3, 3)      :: pp
114      REAL(KIND=dp), DIMENSION(15, 5, 3)      :: pd
115      REAL(KIND=dp), DIMENSION(15, 5, 5)      :: dd
116      ! Derivatives of Rotation Matrices
117      REAL(KIND=dp), DIMENSION(3, 3, 3)   :: sp_d
118      REAL(KIND=dp), DIMENSION(3, 5, 5)   :: sd_d
119      REAL(KIND=dp), DIMENSION(3, 6, 3, 3)   :: pp_d
120      REAL(KIND=dp), DIMENSION(3, 15, 5, 3)   :: pd_d
121      REAL(KIND=dp), DIMENSION(3, 15, 5, 5)   :: dd_d
122   END TYPE rotmat_type
123
124! **************************************************************************************************
125!> \brief  Ewald control type (for periodic SE)
126!> \author Teodoro Laino [tlaino]  - 12.2008
127! **************************************************************************************************
128   TYPE ewald_gks_type
129      REAL(KIND=dp)                            :: alpha
130      TYPE(dg_type), POINTER                   :: dg
131      TYPE(pw_pool_type), POINTER              :: pw_pool
132   END TYPE ewald_gks_type
133
134   TYPE se_int_control_type
135      LOGICAL                                  :: shortrange
136      LOGICAL                                  :: do_ewald_r3
137      LOGICAL                                  :: do_ewald_gks
138      LOGICAL                                  :: pc_coulomb_int
139      INTEGER                                  :: integral_screening
140      INTEGER                                  :: max_multipole
141      TYPE(ewald_gks_type)                     :: ewald_gks
142   END TYPE se_int_control_type
143
144! **************************************************************************************************
145!> \brief Store the value of the tapering function and possibly its derivative
146!>        for screened integrals
147! **************************************************************************************************
148   TYPE se_int_screen_type
149      REAL(KIND=dp)                            :: ft, dft
150   END TYPE se_int_screen_type
151
152! **************************************************************************************************
153!> \brief Taper type use in semi-empirical calculations
154! **************************************************************************************************
155   TYPE se_taper_type
156      TYPE(taper_type), POINTER             :: taper
157      TYPE(taper_type), POINTER             :: taper_cou
158      TYPE(taper_type), POINTER             :: taper_exc
159      TYPE(taper_type), POINTER             :: taper_lrc
160      ! This taper is for KDSO-D integrals
161      TYPE(taper_type), POINTER             :: taper_add
162   END TYPE se_taper_type
163
164   PUBLIC :: semi_empirical_type, &
165             semi_empirical_p_type, &
166             semi_empirical_create, &
167             semi_empirical_release, &
168             rotmat_type, &
169             rotmat_create, &
170             rotmat_release, &
171             get_se_param, &
172             write_se_param, &
173             se_int_control_type, &
174             setup_se_int_control_type, &
175             se_int_screen_type, &
176             se_taper_type, &
177             se_taper_release, &
178             se_taper_create
179
180CONTAINS
181
182! **************************************************************************************************
183!> \brief Allocate semi-empirical type
184!> \param sep ...
185! **************************************************************************************************
186   SUBROUTINE semi_empirical_create(sep)
187      TYPE(semi_empirical_type), POINTER                 :: sep
188
189      CHARACTER(len=*), PARAMETER :: routineN = 'semi_empirical_create', &
190         routineP = moduleN//':'//routineN
191
192      CPASSERT(.NOT. ASSOCIATED(sep))
193      ALLOCATE (sep)
194      ALLOCATE (sep%beta(0:3))
195      ALLOCATE (sep%sto_exponents(0:3))
196      ALLOCATE (sep%zn(0:3))
197      NULLIFY (sep%basis)
198      NULLIFY (sep%w)
199      NULLIFY (sep%w_mpole)
200      NULLIFY (sep%expns3_int)
201      CALL zero_se_param(sep)
202
203   END SUBROUTINE semi_empirical_create
204
205! **************************************************************************************************
206!> \brief Deallocate the semi-empirical type
207!> \param sep ...
208! **************************************************************************************************
209   SUBROUTINE semi_empirical_release(sep)
210
211      TYPE(semi_empirical_type), POINTER                 :: sep
212
213      CHARACTER(len=*), PARAMETER :: routineN = 'semi_empirical_release', &
214         routineP = moduleN//':'//routineN
215
216      INTEGER                                            :: i
217
218      IF (ASSOCIATED(sep)) THEN
219         CALL deallocate_sto_basis_set(sep%basis)
220         CALL semi_empirical_mpole_p_release(sep%w_mpole)
221         IF (ASSOCIATED(sep%beta)) THEN
222            DEALLOCATE (sep%beta)
223         END IF
224         IF (ASSOCIATED(sep%sto_exponents)) THEN
225            DEALLOCATE (sep%sto_exponents)
226         END IF
227         IF (ASSOCIATED(sep%zn)) THEN
228            DEALLOCATE (sep%zn)
229         END IF
230         IF (ASSOCIATED(sep%w)) THEN
231            DEALLOCATE (sep%w)
232         END IF
233         IF (ASSOCIATED(sep%expns3_int)) THEN
234            DO i = 1, SIZE(sep%expns3_int)
235               CALL semi_empirical_expns3_release(sep%expns3_int(i)%expns3)
236            END DO
237            DEALLOCATE (sep%expns3_int)
238         END IF
239         DEALLOCATE (sep)
240      END IF
241
242   END SUBROUTINE semi_empirical_release
243
244! **************************************************************************************************
245!> \brief Zero the whole semi-empirical type
246!> \param sep ...
247! **************************************************************************************************
248   SUBROUTINE zero_se_param(sep)
249      TYPE(semi_empirical_type), POINTER                 :: sep
250
251      CHARACTER(len=*), PARAMETER :: routineN = 'zero_se_param', routineP = moduleN//':'//routineN
252
253      CPASSERT(ASSOCIATED(sep))
254      sep%defined = .FALSE.
255      sep%dorb = .FALSE.
256      sep%extended_basis_set = .FALSE.
257      sep%p_orbitals_on_h = .FALSE.
258      sep%name = ""
259      sep%typ = HUGE(0)
260      sep%core_size = HUGE(0)
261      sep%atm_int_size = HUGE(0)
262      sep%z = HUGE(0)
263      sep%zeff = HUGE(0.0_dp)
264      sep%natorb = 0
265      sep%ngauss = 0
266      sep%eheat = HUGE(0.0_dp)
267
268      sep%zn = 0.0_dp
269      sep%sto_exponents = 0.0_dp
270      sep%beta = 0.0_dp
271
272      sep%uss = 0.0_dp !eV
273      sep%upp = 0.0_dp !eV
274      sep%udd = 0.0_dp !eV
275      sep%uff = 0.0_dp
276      sep%alp = 0.0_dp
277      sep%eisol = 0.0_dp
278      sep%nr = 1
279      sep%acoul = 0.0_dp
280      sep%de = 0.0_dp
281      sep%ass = 0.0_dp
282      sep%asp = 0.0_dp
283      sep%app = 0.0_dp
284      sep%gss = 0.0_dp
285      sep%gsp = 0.0_dp
286      sep%gpp = 0.0_dp
287      sep%gp2 = 0.0_dp
288      sep%gsd = 0.0_dp
289      sep%gpd = 0.0_dp
290      sep%gdd = 0.0_dp
291      sep%hsp = 0.0_dp
292      sep%dd = 0.0_dp
293      sep%qq = 0.0_dp
294      sep%am = 0.0_dp
295      sep%ad = 0.0_dp
296      sep%aq = 0.0_dp
297
298      sep%fn1 = 0.0_dp
299      sep%fn2 = 0.0_dp
300      sep%fn3 = 0.0_dp
301      sep%bfn1 = 0.0_dp
302      sep%bfn2 = 0.0_dp
303      sep%bfn3 = 0.0_dp
304
305      sep%pre = 0.0_dp
306      sep%d = 0.0_dp
307
308      sep%xab = 0.0_dp
309      sep%aab = 0.0_dp
310      sep%a = 0.0_dp
311      sep%b = 0.0_dp
312      sep%c = 0.0_dp
313      sep%rho = 0.0_dp
314
315      sep%f0dd = 0.0_dp
316      sep%f2dd = 0.0_dp
317      sep%f4dd = 0.0_dp
318      sep%f0sd = 0.0_dp
319      sep%f0pd = 0.0_dp
320      sep%f2pd = 0.0_dp
321      sep%g1pd = 0.0_dp
322      sep%g2sd = 0.0_dp
323      sep%g3pd = 0.0_dp
324      sep%ko = 0.0_dp
325      sep%cs = 0.0_dp
326      sep%onec2el = 0.0_dp
327
328   END SUBROUTINE zero_se_param
329
330! **************************************************************************************************
331!> \brief Get info from the semi-empirical type
332!> \param sep ...
333!> \param name ...
334!> \param typ ...
335!> \param defined ...
336!> \param z ...
337!> \param zeff ...
338!> \param natorb ...
339!> \param eheat ...
340!> \param beta ...
341!> \param sto_exponents ...
342!> \param uss ...
343!> \param upp ...
344!> \param udd ...
345!> \param uff ...
346!> \param alp ...
347!> \param eisol ...
348!> \param gss ...
349!> \param gsp ...
350!> \param gpp ...
351!> \param gp2 ...
352!> \param acoul ...
353!> \param nr ...
354!> \param de ...
355!> \param ass ...
356!> \param asp ...
357!> \param app ...
358!> \param hsp ...
359!> \param gsd ...
360!> \param gpd ...
361!> \param gdd ...
362!> \param ppddg ...
363!> \param dpddg ...
364!> \param ngauss ...
365! **************************************************************************************************
366   SUBROUTINE get_se_param(sep, name, typ, defined, z, zeff, natorb, eheat, &
367                           beta, sto_exponents, uss, upp, udd, uff, alp, eisol, gss, gsp, gpp, gp2, &
368                           acoul, nr, de, ass, asp, app, hsp, gsd, gpd, gdd, ppddg, dpddg, ngauss)
369
370      TYPE(semi_empirical_type), POINTER                 :: sep
371      CHARACTER(LEN=default_string_length), &
372         INTENT(OUT), OPTIONAL                           :: name
373      INTEGER, INTENT(OUT), OPTIONAL                     :: typ
374      LOGICAL, INTENT(OUT), OPTIONAL                     :: defined
375      INTEGER, INTENT(OUT), OPTIONAL                     :: z
376      REAL(KIND=dp), INTENT(OUT), OPTIONAL               :: zeff
377      INTEGER, INTENT(OUT), OPTIONAL                     :: natorb
378      REAL(KIND=dp), OPTIONAL                            :: eheat
379      REAL(KIND=dp), DIMENSION(:), OPTIONAL, POINTER     :: beta, sto_exponents
380      REAL(KIND=dp), OPTIONAL                            :: uss, upp, udd, uff, alp, eisol, gss, &
381                                                            gsp, gpp, gp2, acoul
382      INTEGER, INTENT(OUT), OPTIONAL                     :: nr
383      REAL(KIND=dp), OPTIONAL                            :: de, ass, asp, app, hsp, gsd, gpd, gdd
384      REAL(KIND=dp), DIMENSION(2), OPTIONAL              :: ppddg, dpddg
385      INTEGER, INTENT(OUT), OPTIONAL                     :: ngauss
386
387      CHARACTER(len=*), PARAMETER :: routineN = 'get_se_param', routineP = moduleN//':'//routineN
388
389      IF (ASSOCIATED(sep)) THEN
390         IF (PRESENT(name)) name = sep%name
391         IF (PRESENT(typ)) typ = sep%typ
392         IF (PRESENT(defined)) defined = sep%defined
393         IF (PRESENT(z)) z = sep%z
394         IF (PRESENT(zeff)) zeff = sep%zeff
395         IF (PRESENT(natorb)) natorb = sep%natorb
396         IF (PRESENT(eheat)) eheat = sep%eheat
397         IF (PRESENT(beta)) beta => sep%beta
398         IF (PRESENT(sto_exponents)) sto_exponents => sep%sto_exponents
399         IF (PRESENT(ngauss)) ngauss = sep%ngauss
400         IF (PRESENT(uss)) uss = sep%uss
401         IF (PRESENT(upp)) upp = sep%upp
402         IF (PRESENT(udd)) udd = sep%udd
403         IF (PRESENT(uff)) uff = sep%uff
404         IF (PRESENT(alp)) alp = sep%alp
405         IF (PRESENT(eisol)) eisol = sep%eisol
406         IF (PRESENT(nr)) nr = sep%nr
407         IF (PRESENT(acoul)) acoul = sep%acoul
408         IF (PRESENT(de)) de = sep%de
409         IF (PRESENT(ass)) ass = sep%ass
410         IF (PRESENT(asp)) asp = sep%asp
411         IF (PRESENT(app)) app = sep%app
412         IF (PRESENT(gss)) gss = sep%gss
413         IF (PRESENT(gsp)) gsp = sep%gsp
414         IF (PRESENT(gpp)) gpp = sep%gpp
415         IF (PRESENT(gp2)) gp2 = sep%gp2
416         IF (PRESENT(hsp)) hsp = sep%hsp
417         IF (PRESENT(gsd)) gsd = sep%gsd
418         IF (PRESENT(gpd)) gpd = sep%gpd
419         IF (PRESENT(gdd)) gdd = sep%gdd
420         IF (PRESENT(ppddg)) ppddg = sep%pre
421         IF (PRESENT(dpddg)) dpddg = sep%d
422      ELSE
423         CPABORT("The pointer sep is not associated")
424      END IF
425
426   END SUBROUTINE get_se_param
427
428! **************************************************************************************************
429!> \brief Set info from the semi-empirical type
430!> \param sep ...
431!> \param name ...
432!> \param typ ...
433!> \param defined ...
434!> \param z ...
435!> \param zeff ...
436!> \param natorb ...
437!> \param eheat ...
438!> \param beta ...
439!> \param sto_exponents ...
440!> \param uss ...
441!> \param upp ...
442!> \param udd ...
443!> \param uff ...
444!> \param alp ...
445!> \param eisol ...
446!> \param gss ...
447!> \param gsp ...
448!> \param gpp ...
449!> \param gp2 ...
450!> \param acoul ...
451!> \param nr ...
452!> \param de ...
453!> \param ass ...
454!> \param asp ...
455!> \param app ...
456!> \param hsp ...
457!> \param gsd ...
458!> \param gpd ...
459!> \param gdd ...
460!> \param ppddg ...
461!> \param dpddg ...
462!> \param ngauss ...
463! **************************************************************************************************
464   SUBROUTINE set_se_param(sep, name, typ, defined, z, zeff, natorb, eheat, &
465                           beta, sto_exponents, uss, upp, udd, uff, alp, eisol, gss, gsp, gpp, gp2, &
466                           acoul, nr, de, ass, asp, app, hsp, gsd, gpd, gdd, ppddg, dpddg, ngauss)
467
468      TYPE(semi_empirical_type), POINTER                 :: sep
469      CHARACTER(LEN=default_string_length), INTENT(IN), &
470         OPTIONAL                                        :: name
471      INTEGER, INTENT(IN), OPTIONAL                      :: typ
472      LOGICAL, INTENT(IN), OPTIONAL                      :: defined
473      INTEGER, INTENT(IN), OPTIONAL                      :: z
474      REAL(KIND=dp), INTENT(IN), OPTIONAL                :: zeff
475      INTEGER, INTENT(IN), OPTIONAL                      :: natorb
476      REAL(KIND=dp), OPTIONAL                            :: eheat
477      REAL(dp), DIMENSION(0:), OPTIONAL                  :: beta
478      REAL(KIND=dp), DIMENSION(:), OPTIONAL              :: sto_exponents
479      REAL(KIND=dp), OPTIONAL                            :: uss, upp, udd, uff, alp, eisol, gss, &
480                                                            gsp, gpp, gp2, acoul
481      INTEGER, INTENT(IN), OPTIONAL                      :: nr
482      REAL(KIND=dp), OPTIONAL                            :: de, ass, asp, app, hsp, gsd, gpd, gdd
483      REAL(dp), DIMENSION(2), OPTIONAL                   :: ppddg, dpddg
484      INTEGER, INTENT(IN), OPTIONAL                      :: ngauss
485
486      CHARACTER(len=*), PARAMETER :: routineN = 'set_se_param', routineP = moduleN//':'//routineN
487
488      IF (ASSOCIATED(sep)) THEN
489         IF (PRESENT(name)) sep%name = name
490         IF (PRESENT(typ)) sep%typ = typ
491         IF (PRESENT(defined)) sep%defined = defined
492         IF (PRESENT(z)) sep%z = z
493         IF (PRESENT(zeff)) sep%zeff = zeff
494         IF (PRESENT(natorb)) sep%natorb = natorb
495         IF (PRESENT(eheat)) sep%eheat = eheat
496         IF (PRESENT(beta)) sep%beta = beta
497         IF (PRESENT(sto_exponents)) sep%sto_exponents = sto_exponents
498         IF (PRESENT(ngauss)) sep%ngauss = ngauss
499         IF (PRESENT(uss)) sep%uss = uss
500         IF (PRESENT(upp)) sep%upp = upp
501         IF (PRESENT(udd)) sep%udd = udd
502         IF (PRESENT(uff)) sep%uff = uff
503         IF (PRESENT(alp)) sep%alp = alp
504         IF (PRESENT(eisol)) sep%eisol = eisol
505         IF (PRESENT(acoul)) sep%acoul = acoul
506         IF (PRESENT(nr)) sep%nr = nr
507         IF (PRESENT(de)) sep%de = de
508         IF (PRESENT(ass)) sep%ass = ass
509         IF (PRESENT(asp)) sep%asp = asp
510         IF (PRESENT(app)) sep%app = app
511         IF (PRESENT(gss)) sep%gss = gss
512         IF (PRESENT(gsp)) sep%gsp = gsp
513         IF (PRESENT(gpp)) sep%gpp = gpp
514         IF (PRESENT(gp2)) sep%gp2 = gp2
515         IF (PRESENT(hsp)) sep%hsp = hsp
516         IF (PRESENT(gsd)) sep%gsd = gsd
517         IF (PRESENT(gpd)) sep%gpd = gpd
518         IF (PRESENT(gdd)) sep%gdd = gdd
519         IF (PRESENT(ppddg)) sep%pre = ppddg
520         IF (PRESENT(dpddg)) sep%d = dpddg
521      ELSE
522         CPABORT("The pointer sep is not associated")
523      END IF
524
525   END SUBROUTINE set_se_param
526
527! **************************************************************************************************
528!> \brief Creates rotmat type
529!> \param rotmat ...
530! **************************************************************************************************
531   SUBROUTINE rotmat_create(rotmat)
532      TYPE(rotmat_type), POINTER                         :: rotmat
533
534      CHARACTER(len=*), PARAMETER :: routineN = 'rotmat_create', routineP = moduleN//':'//routineN
535
536      CPASSERT(.NOT. ASSOCIATED(rotmat))
537      ALLOCATE (rotmat)
538
539   END SUBROUTINE rotmat_create
540
541! **************************************************************************************************
542!> \brief Releases rotmat type
543!> \param rotmat ...
544! **************************************************************************************************
545   SUBROUTINE rotmat_release(rotmat)
546      TYPE(rotmat_type), POINTER                         :: rotmat
547
548      CHARACTER(len=*), PARAMETER :: routineN = 'rotmat_release', routineP = moduleN//':'//routineN
549
550      IF (ASSOCIATED(rotmat)) THEN
551         DEALLOCATE (rotmat)
552      END IF
553
554   END SUBROUTINE rotmat_release
555
556! **************************************************************************************************
557!> \brief Setup the Semiempirical integral control type
558!> \param se_int_control ...
559!> \param shortrange ...
560!> \param do_ewald_r3 ...
561!> \param do_ewald_gks ...
562!> \param integral_screening ...
563!> \param max_multipole ...
564!> \param pc_coulomb_int ...
565!> \author Teodoro Laino [tlaino] - 12.2008
566! **************************************************************************************************
567   SUBROUTINE setup_se_int_control_type(se_int_control, shortrange, do_ewald_r3, &
568                                        do_ewald_gks, integral_screening, max_multipole, pc_coulomb_int)
569      TYPE(se_int_control_type)                          :: se_int_control
570      LOGICAL, INTENT(IN)                                :: shortrange, do_ewald_r3, do_ewald_gks
571      INTEGER, INTENT(IN)                                :: integral_screening, max_multipole
572      LOGICAL, INTENT(IN)                                :: pc_coulomb_int
573
574      se_int_control%shortrange = shortrange
575      se_int_control%do_ewald_r3 = do_ewald_r3
576      se_int_control%integral_screening = integral_screening
577      ! This makes the assignment independent of the value of the different constants
578      SELECT CASE (max_multipole)
579      CASE (do_multipole_none)
580         se_int_control%max_multipole = -1
581      CASE (do_multipole_charge)
582         se_int_control%max_multipole = 0
583      CASE (do_multipole_dipole)
584         se_int_control%max_multipole = 1
585      CASE (do_multipole_quadrupole)
586         se_int_control%max_multipole = 2
587      END SELECT
588
589      se_int_control%do_ewald_gks = do_ewald_gks
590      se_int_control%pc_coulomb_int = pc_coulomb_int
591      NULLIFY (se_int_control%ewald_gks%dg, se_int_control%ewald_gks%pw_pool)
592
593   END SUBROUTINE setup_se_int_control_type
594
595! **************************************************************************************************
596!> \brief Creates the taper type used in SE calculations
597!> \param se_taper ...
598!> \param integral_screening ...
599!> \param do_ewald ...
600!> \param taper_cou ...
601!> \param range_cou ...
602!> \param taper_exc ...
603!> \param range_exc ...
604!> \param taper_scr ...
605!> \param range_scr ...
606!> \param taper_lrc ...
607!> \param range_lrc ...
608!> \author Teodoro Laino [tlaino] - 03.2009
609! **************************************************************************************************
610   SUBROUTINE se_taper_create(se_taper, integral_screening, do_ewald, &
611                              taper_cou, range_cou, taper_exc, range_exc, taper_scr, range_scr, &
612                              taper_lrc, range_lrc)
613      TYPE(se_taper_type), POINTER                       :: se_taper
614      INTEGER, INTENT(IN)                                :: integral_screening
615      LOGICAL, INTENT(IN)                                :: do_ewald
616      REAL(KIND=dp), INTENT(IN)                          :: taper_cou, range_cou, taper_exc, &
617                                                            range_exc, taper_scr, range_scr, &
618                                                            taper_lrc, range_lrc
619
620      CHARACTER(len=*), PARAMETER :: routineN = 'se_taper_create', &
621         routineP = moduleN//':'//routineN
622
623      CPASSERT(.NOT. ASSOCIATED(se_taper))
624      ALLOCATE (se_taper)
625      NULLIFY (se_taper%taper)
626      NULLIFY (se_taper%taper_cou)
627      NULLIFY (se_taper%taper_exc)
628      NULLIFY (se_taper%taper_lrc)
629      NULLIFY (se_taper%taper_add)
630      ! Create the sub-typo taper
631      CALL taper_create(se_taper%taper_cou, taper_cou, range_cou)
632      CALL taper_create(se_taper%taper_exc, taper_exc, range_exc)
633      IF (integral_screening == do_se_IS_kdso_d) THEN
634         CALL taper_create(se_taper%taper_add, taper_scr, range_scr)
635      END IF
636      IF ((integral_screening /= do_se_IS_slater) .AND. do_ewald) THEN
637         CALL taper_create(se_taper%taper_lrc, taper_lrc, range_lrc)
638      END IF
639   END SUBROUTINE se_taper_create
640
641! **************************************************************************************************
642!> \brief Releases the taper type used in SE calculations
643!> \param se_taper ...
644!> \author Teodoro Laino [tlaino] - 03.2009
645! **************************************************************************************************
646   SUBROUTINE se_taper_release(se_taper)
647      TYPE(se_taper_type), POINTER                       :: se_taper
648
649      CHARACTER(len=*), PARAMETER :: routineN = 'se_taper_release', &
650         routineP = moduleN//':'//routineN
651
652      IF (ASSOCIATED(se_taper)) THEN
653         CALL taper_release(se_taper%taper_cou)
654         CALL taper_release(se_taper%taper_exc)
655         CALL taper_release(se_taper%taper_lrc)
656         CALL taper_release(se_taper%taper_add)
657
658         DEALLOCATE (se_taper)
659      END IF
660   END SUBROUTINE se_taper_release
661
662! **************************************************************************************************
663!> \brief Writes the semi-empirical type
664!> \param sep ...
665!> \param subsys_section ...
666!> \par History
667!>        04.2008 Teodoro Laino [tlaino] - University of Zurich: rewriting with
668!>                support for the whole set of parameters
669! **************************************************************************************************
670   SUBROUTINE write_se_param(sep, subsys_section)
671
672      TYPE(semi_empirical_type), POINTER                 :: sep
673      TYPE(section_vals_type), POINTER                   :: subsys_section
674
675      CHARACTER(len=*), PARAMETER :: routineN = 'write_se_param', routineP = moduleN//':'//routineN
676      CHARACTER(LEN=1), DIMENSION(0:3), PARAMETER :: orb_lab = (/"S", "P", "D", "F"/)
677      CHARACTER(LEN=2), DIMENSION(0:3), PARAMETER :: z_lab = (/"ZS", "ZP", "ZD", "ZF"/)
678      CHARACTER(LEN=3), DIMENSION(0:3), PARAMETER :: zeta_lab = (/"ZSN", "ZPN", "ZDN", "ZFN"/)
679      CHARACTER(LEN=5), DIMENSION(0:3), PARAMETER :: &
680         beta_lab = (/"BETAS", "BETAP", "BETAD", "BETAF"/)
681      CHARACTER(LEN=default_string_length)               :: i_string, name
682      INTEGER                                            :: i, l, natorb, ngauss, nr, output_unit, &
683                                                            typ, z
684      LOGICAL                                            :: defined
685      REAL(KIND=dp)                                      :: acoul, alp, app, asp, ass, de, eheat, &
686                                                            eisol, gp2, gpp, gsp, gss, hsp, udd, &
687                                                            uff, upp, uss, zeff
688      CHARACTER(LEN=3), DIMENSION(0:3), PARAMETER :: u_lab = (/"USS", "UPP", "UDD", "UFF"/)
689
690      REAL(KIND=dp), DIMENSION(0:3)                      :: u
691      REAL(KIND=dp), DIMENSION(2)                        :: dpddg, ppddg
692      REAL(KIND=dp), DIMENSION(:), POINTER               :: beta, sexp
693      TYPE(cp_logger_type), POINTER                      :: logger
694
695      NULLIFY (logger)
696      logger => cp_get_default_logger()
697      IF (ASSOCIATED(sep) .AND. BTEST(cp_print_key_should_output(logger%iter_info, subsys_section, &
698                                                                 "PRINT%KINDS/SE_PARAMETERS"), cp_p_file)) THEN
699
700         output_unit = cp_print_key_unit_nr(logger, subsys_section, "PRINT%KINDS/SE_PARAMETERS", &
701                                            extension=".Log")
702
703         IF (output_unit > 0) THEN
704            CALL get_se_param(sep, name=name, typ=typ, defined=defined, &
705                              z=z, zeff=zeff, natorb=natorb, eheat=eheat, beta=beta, &
706                              sto_exponents=sexp, uss=uss, upp=upp, udd=udd, uff=uff, &
707                              alp=alp, eisol=eisol, gss=gss, gsp=gsp, gpp=gpp, gp2=gp2, &
708                              de=de, ass=ass, asp=asp, app=app, hsp=hsp, ppddg=ppddg, &
709                              acoul=acoul, nr=nr, dpddg=dpddg, ngauss=ngauss)
710
711            u(0) = uss
712            u(1) = upp
713            u(2) = udd
714            u(3) = uff
715
716            SELECT CASE (typ)
717            CASE DEFAULT
718               CPABORT("Semiempirical method unknown")
719            CASE (do_method_am1)
720               WRITE (UNIT=output_unit, FMT="(/,A,T35,A,T67,A14)") &
721                  " Semi empirical parameters: ", "Austin Model 1 (AM1)", TRIM(name)
722            CASE (do_method_rm1)
723               WRITE (UNIT=output_unit, FMT="(/,A,T35,A,T67,A14)") &
724                  " Semi empirical parameters: ", "Recife Model 1 (RM1)", TRIM(name)
725            CASE (do_method_pm3)
726               WRITE (UNIT=output_unit, FMT="(/,A,T35,A,T67,A14)") &
727                  " Semi empirical parameters: ", "Parametric Method 3 (PM3) ", TRIM(name)
728            CASE (do_method_pnnl)
729               WRITE (UNIT=output_unit, FMT="(/,A,T35,A,T67,A14)") &
730                  " Semi empirical parameters: ", "PNNL method ", TRIM(name)
731            CASE (do_method_pm6)
732               WRITE (UNIT=output_unit, FMT="(/,A,T35,A,T67,A14)") &
733                  " Semi empirical parameters: ", "Parametric Method 6 (PM6) ", TRIM(name)
734            CASE (do_method_pm6fm)
735               WRITE (UNIT=output_unit, FMT="(/,A,T35,A,T67,A14)") &
736                  " Semi empirical parameters: ", "Parametric Method 6 (PM6-FM) ", TRIM(name)
737            CASE (do_method_pdg)
738               WRITE (UNIT=output_unit, FMT="(/,A,T35,A,T67,A14)") &
739                  " Semi empirical parameters: ", "PDDG/PM3 ", TRIM(name)
740            CASE (do_method_mndo)
741               WRITE (UNIT=output_unit, FMT="(/,A,T35,A,T67,A14)") &
742                  " Semi empirical parameters: ", "MNDO ", TRIM(name)
743            CASE (do_method_mndod)
744               WRITE (UNIT=output_unit, FMT="(/,A,T35,A,T67,A14)") &
745                  " Semi empirical parameters: ", "MNDOD", TRIM(name)
746            END SELECT
747
748            ! If defined print all its semi-empirical parameters
749            IF (defined) THEN
750               WRITE (UNIT=output_unit, FMT="(T16,A,T71,F10.2)") &
751                  "Effective core charge:", zeff
752               WRITE (UNIT=output_unit, FMT="(T16,A,T71,I10)") &
753                  "Number of orbitals:", natorb, &
754                  "Basis set expansion (STO-NG)", ngauss
755               WRITE (UNIT=output_unit, FMT="(T16,A,T66,F15.5)") &
756                  "Atomic heat of formation [kcal/mol]:", eheat*kcalmol
757               DO l = 0, 3
758                  IF (ABS(beta(l)) > 0._dp) THEN
759                     WRITE (UNIT=output_unit, FMT="(T16,A,I2)") "Parameters for Shell: ", l
760                     WRITE (UNIT=output_unit, FMT="(T22,A5,T30,A,T64,F17.4)") &
761                        ADJUSTR(z_lab(l)), "- "//"Slater  Exponent for "//orb_lab(l)//"     [A]: ", sexp(l)
762                     WRITE (UNIT=output_unit, FMT="(T22,A5,T30,A,T64,F17.4)") &
763                        ADJUSTR(u_lab(l)), "- "//"One Center Energy for "//orb_lab(l)//"   [eV]: ", u(l)*evolt
764                     WRITE (UNIT=output_unit, FMT="(T22,A5,T30,A,T64,F17.4)") &
765                        ADJUSTR(beta_lab(l)), "- "//"Beta Parameter for "//orb_lab(l)//"      [eV]: ", beta(l)*evolt
766                     WRITE (UNIT=output_unit, FMT="(T22,A5,T30,A,T64,F17.4)") &
767                        ADJUSTR(zeta_lab(l)), "- "//"Internal Exponent for "//orb_lab(l)//" [a.u.]: ", sep%zn(l)
768                  END IF
769               END DO
770               WRITE (UNIT=output_unit, FMT="(/,T16,A)") "Additional Parameters (Derived or Fitted):"
771               WRITE (UNIT=output_unit, FMT="(T16,A11,T30,A,T69,F12.4)") &
772                  ADJUSTR("ALP"), "- "//"Alpha Parameter for Core    [A^-1]: ", alp/angstrom
773               WRITE (UNIT=output_unit, FMT="(T16,A11,T30,A,T69,F12.4)") &
774                  ADJUSTR("EISOL"), "- "//"Atomic Energy (Calculated)    [eV]: ", eisol*evolt
775               ! One center Two electron Integrals
776               WRITE (UNIT=output_unit, FMT="(T16,A11,T30,A,T69,F12.4)") &
777                  ADJUSTR("GSS"), "- "//"One Center Integral (SS ,SS ) [eV]: ", gss*evolt
778               WRITE (UNIT=output_unit, FMT="(T16,A11,T30,A,T69,F12.4)") &
779                  ADJUSTR("GSP"), "- "//"One Center Integral (SS ,PP ) [eV]: ", gsp*evolt
780               WRITE (UNIT=output_unit, FMT="(T16,A11,T30,A,T69,F12.4)") &
781                  ADJUSTR("GPP"), "- "//"One Center Integral (PP ,PP ) [eV]: ", gpp*evolt
782               WRITE (UNIT=output_unit, FMT="(T16,A11,T30,A,T69,F12.4)") &
783                  ADJUSTR("GP2"), "- "//"One Center Integral (PP*,PP*) [eV]: ", gp2*evolt
784               WRITE (UNIT=output_unit, FMT="(T16,A11,T30,A,T69,F12.4)") &
785                  ADJUSTR("HSP"), "- "//"One Center Integral (SP ,SP ) [eV]: ", hsp*evolt
786               ! Slater Condon Parameters
787               IF (sep%dorb) THEN
788                  WRITE (UNIT=output_unit, FMT="(T16,A11,T30,A,T69,F12.4)") &
789                     ADJUSTR("F0DD"), "- "//"Slater Condon Parameter F0DD  [eV]: ", sep%f0dd
790                  WRITE (UNIT=output_unit, FMT="(T16,A11,T30,A,T69,F12.4)") &
791                     ADJUSTR("F2DD"), "- "//"Slater Condon Parameter F2DD  [eV]: ", sep%f2dd
792                  WRITE (UNIT=output_unit, FMT="(T16,A11,T30,A,T69,F12.4)") &
793                     ADJUSTR("F4DD"), "- "//"Slater Condon Parameter F4DD  [eV]: ", sep%f4dd
794                  WRITE (UNIT=output_unit, FMT="(T16,A11,T30,A,T69,F12.4)") &
795                     ADJUSTR("FOSD"), "- "//"Slater Condon Parameter FOSD  [eV]: ", sep%f0sd
796                  WRITE (UNIT=output_unit, FMT="(T16,A11,T30,A,T69,F12.4)") &
797                     ADJUSTR("G2SD"), "- "//"Slater Condon Parameter G2SD  [eV]: ", sep%g2sd
798                  WRITE (UNIT=output_unit, FMT="(T16,A11,T30,A,T69,F12.4)") &
799                     ADJUSTR("F0PD"), "- "//"Slater Condon Parameter F0PD  [eV]: ", sep%f0pd
800                  WRITE (UNIT=output_unit, FMT="(T16,A11,T30,A,T69,F12.4)") &
801                     ADJUSTR("F2PD"), "- "//"Slater Condon Parameter F2PD  [eV]: ", sep%f2pd
802                  WRITE (UNIT=output_unit, FMT="(T16,A11,T30,A,T69,F12.4)") &
803                     ADJUSTR("G1PD"), "- "//"Slater Condon Parameter G1PD  [eV]: ", sep%g1pd
804                  WRITE (UNIT=output_unit, FMT="(T16,A11,T30,A,T69,F12.4)") &
805                     ADJUSTR("G3PD"), "- "//"Slater Condon Parameter G3PD  [eV]: ", sep%g3pd
806               END IF
807               ! Charge Separation
808               WRITE (UNIT=output_unit, FMT="(T16,A11,T30,A,T69,F12.4)") &
809                  ADJUSTR("DD2"), "- "//"Charge Separation  SP, L=1  [bohr]: ", sep%cs(2)
810               WRITE (UNIT=output_unit, FMT="(T16,A11,T30,A,T69,F12.4)") &
811                  ADJUSTR("DD3"), "- "//"Charge Separation  PP, L=2  [bohr]: ", sep%cs(3)
812               IF (sep%dorb) THEN
813                  WRITE (UNIT=output_unit, FMT="(T16,A11,T30,A,T69,F12.4)") &
814                     ADJUSTR("DD4"), "- "//"Charge Separation  SD, L=2  [bohr]: ", sep%cs(4)
815                  WRITE (UNIT=output_unit, FMT="(T16,A11,T30,A,T69,F12.4)") &
816                     ADJUSTR("DD5"), "- "//"Charge Separation  PD, L=1  [bohr]: ", sep%cs(5)
817                  WRITE (UNIT=output_unit, FMT="(T16,A11,T30,A,T69,F12.4)") &
818                     ADJUSTR("DD6"), "- "//"Charge Separation  DD, L=2  [bohr]: ", sep%cs(6)
819               END IF
820               ! Klopman-Ohno Terms
821               WRITE (UNIT=output_unit, FMT="(T16,A11,T30,A,T69,F12.4)") &
822                  ADJUSTR("PO1"), "- "//"Klopman-Ohno term, SS, L=0  [bohr]: ", sep%ko(1)
823               WRITE (UNIT=output_unit, FMT="(T16,A11,T30,A,T69,F12.4)") &
824                  ADJUSTR("PO2"), "- "//"Klopman-Ohno term, SP, L=1  [bohr]: ", sep%ko(2)
825               WRITE (UNIT=output_unit, FMT="(T16,A11,T30,A,T69,F12.4)") &
826                  ADJUSTR("PO3"), "- "//"Klopman-Ohno term, PP, L=2  [bohr]: ", sep%ko(3)
827               IF (sep%dorb) THEN
828                  WRITE (UNIT=output_unit, FMT="(T16,A11,T30,A,T69,F12.4)") &
829                     ADJUSTR("PO4"), "- "//"Klopman-Ohno term, SD, L=2  [bohr]: ", sep%ko(4)
830                  WRITE (UNIT=output_unit, FMT="(T16,A11,T30,A,T69,F12.4)") &
831                     ADJUSTR("PO5"), "- "//"Klopman-Ohno term, PD, L=1  [bohr]: ", sep%ko(5)
832                  WRITE (UNIT=output_unit, FMT="(T16,A11,T30,A,T69,F12.4)") &
833                     ADJUSTR("PO6"), "- "//"Klopman-Ohno term, DD, L=2  [bohr]: ", sep%ko(6)
834                  WRITE (UNIT=output_unit, FMT="(T16,A11,T30,A,T69,F12.4)") &
835                     ADJUSTR("PO7"), "- "//"Klopman-Ohno term, PP, L=0  [bohr]: ", sep%ko(7)
836                  WRITE (UNIT=output_unit, FMT="(T16,A11,T30,A,T69,F12.4)") &
837                     ADJUSTR("PO8"), "- "//"Klopman-Ohno term, DD, L=0  [bohr]: ", sep%ko(8)
838               END IF
839               WRITE (UNIT=output_unit, FMT="(T16,A11,T30,A,T69,F12.4)") &
840                  ADJUSTR("PO9"), "- "//"Klopman-Ohno term, CORE     [bohr]: ", sep%ko(9)
841               SELECT CASE (typ)
842               CASE (do_method_am1, do_method_rm1, do_method_pm3, do_method_pdg, do_method_pnnl)
843                  IF (typ == do_method_pnnl) THEN
844                     WRITE (UNIT=output_unit, FMT="(T16,A11,T30,A,T69,F12.4)") &
845                        ADJUSTR("ASS"), "- "//" SS polarization [au]: ", sep%ass
846                     WRITE (UNIT=output_unit, FMT="(T16,A11,T30,A,T69,F12.4)") &
847                        ADJUSTR("ASP"), "- "//" SP polarization [au]: ", sep%asp
848                     WRITE (UNIT=output_unit, FMT="(T16,A11,T30,A,T69,F12.4)") &
849                        ADJUSTR("APP"), "- "//" PP polarization[au]: ", sep%app
850                     WRITE (UNIT=output_unit, FMT="(T16,A11,T30,A,T69,F12.4)") &
851                        ADJUSTR("DE"), "- "//" Dispersion Parameter [eV]: ", sep%de*evolt
852                     WRITE (UNIT=output_unit, FMT="(T16,A11,T30,A,T69,F12.4)") &
853                        ADJUSTR("ACOUL"), "- "//" Slater parameter: ", sep%acoul
854                     WRITE (UNIT=output_unit, FMT="(T16,A11,T30,A,T69,I12)") &
855                        ADJUSTR("NR"), "- "//" Slater parameter: ", sep%nr
856                  ELSEIF ((typ == do_method_am1 .OR. typ == do_method_rm1) .AND. sep%z == 5) THEN
857                     ! Standard case
858                     DO i = 1, SIZE(sep%bfn1, 1)
859                        i_string = cp_to_string(i)
860                        WRITE (UNIT=output_unit, FMT="(T16,A11,T30,A,T69,F12.4)") &
861                           ADJUSTR("FN1"//TRIM(ADJUSTL(i_string))//"_ALL"), &
862                           "- "//"Core-Core VDW, Multiplier   [a.u.]: ", sep%bfn1(i, 1)
863                        WRITE (UNIT=output_unit, FMT="(T16,A11,T30,A,T69,F12.4)") &
864                           ADJUSTR("FN2"//TRIM(ADJUSTL(i_string))//"_ALL"), &
865                           "- "//"Core-Core VDW, Exponent     [a.u.]: ", sep%bfn2(i, 1)
866                        WRITE (UNIT=output_unit, FMT="(T16,A11,T30,A,T69,F12.4)") &
867                           ADJUSTR("FN3"//TRIM(ADJUSTL(i_string))//"_ALL"), &
868                           "- "//"Core-Core VDW, Position     [a.u.]: ", sep%bfn3(i, 1)
869                     END DO
870                     ! Special Case : Hydrogen
871                     DO i = 1, SIZE(sep%bfn1, 1)
872                        i_string = cp_to_string(i)
873                        WRITE (UNIT=output_unit, FMT="(T16,A11,T30,A,T69,F12.4)") &
874                           ADJUSTR("FN1"//TRIM(ADJUSTL(i_string))//"_H"), &
875                           "- "//"Core-Core VDW, Multiplier   [a.u.]: ", sep%bfn1(i, 2)
876                        WRITE (UNIT=output_unit, FMT="(T16,A11,T30,A,T69,F12.4)") &
877                           ADJUSTR("FN2"//TRIM(ADJUSTL(i_string))//"_H"), &
878                           "- "//"Core-Core VDW, Exponent     [a.u.]: ", sep%bfn2(i, 2)
879                        WRITE (UNIT=output_unit, FMT="(T16,A11,T30,A,T69,F12.4)") &
880                           ADJUSTR("FN3"//TRIM(ADJUSTL(i_string))//"_H"), &
881                           "- "//"Core-Core VDW, Position     [a.u.]: ", sep%bfn3(i, 2)
882                     END DO
883                     ! Special Case : Carbon
884                     DO i = 1, SIZE(sep%bfn1, 1)
885                        i_string = cp_to_string(i)
886                        WRITE (UNIT=output_unit, FMT="(T16,A11,T30,A,T69,F12.4)") &
887                           ADJUSTR("FN1"//TRIM(ADJUSTL(i_string))//"_C"), &
888                           "- "//"Core-Core VDW, Multiplier   [a.u.]: ", sep%bfn1(i, 3)
889                        WRITE (UNIT=output_unit, FMT="(T16,A11,T30,A,T69,F12.4)") &
890                           ADJUSTR("FN2"//TRIM(ADJUSTL(i_string))//"_C"), &
891                           "- "//"Core-Core VDW, Exponent     [a.u.]: ", sep%bfn2(i, 3)
892                        WRITE (UNIT=output_unit, FMT="(T16,A11,T30,A,T69,F12.4)") &
893                           ADJUSTR("FN3"//TRIM(ADJUSTL(i_string))//"_C"), &
894                           "- "//"Core-Core VDW, Position     [a.u.]: ", sep%bfn3(i, 3)
895                     END DO
896                     ! Special Case : Halogens
897                     DO i = 1, SIZE(sep%bfn1, 1)
898                        i_string = cp_to_string(i)
899                        WRITE (UNIT=output_unit, FMT="(T16,A11,T30,A,T69,F12.4)") &
900                           ADJUSTR("FN1"//TRIM(ADJUSTL(i_string))//"_HALO"), &
901                           "- "//"Core-Core VDW, Multiplier   [a.u.]: ", sep%bfn1(i, 4)
902                        WRITE (UNIT=output_unit, FMT="(T16,A11,T30,A,T69,F12.4)") &
903                           ADJUSTR("FN2"//TRIM(ADJUSTL(i_string))//"_HALO"), &
904                           "- "//"Core-Core VDW, Exponent     [a.u.]: ", sep%bfn2(i, 4)
905                        WRITE (UNIT=output_unit, FMT="(T16,A11,T30,A,T69,F12.4)") &
906                           ADJUSTR("FN3"//TRIM(ADJUSTL(i_string))//"_HALO"), &
907                           "- "//"Core-Core VDW, Position     [a.u.]: ", sep%bfn3(i, 4)
908                     END DO
909                  ELSE
910                     DO i = 1, SIZE(sep%fn1, 1)
911                        i_string = cp_to_string(i)
912                        ! Skip the printing of params that are zero..
913                        IF (sep%fn1(i) == 0.0_dp .AND. sep%fn2(i) == 0.0_dp .AND. sep%fn3(i) == 0.0_dp) CYCLE
914                        WRITE (UNIT=output_unit, FMT="(T16,A11,T30,A,T69,F12.4)") &
915                           ADJUSTR("FN1"//TRIM(ADJUSTL(i_string))), &
916                           "- "//"Core-Core VDW, Multiplier   [a.u.]: ", sep%fn1(i)
917                        WRITE (UNIT=output_unit, FMT="(T16,A11,T30,A,T69,F12.4)") &
918                           ADJUSTR("FN2"//TRIM(ADJUSTL(i_string))), &
919                           "- "//"Core-Core VDW, Exponent     [a.u.]: ", sep%fn2(i)
920                        WRITE (UNIT=output_unit, FMT="(T16,A11,T30,A,T69,F12.4)") &
921                           ADJUSTR("FN3"//TRIM(ADJUSTL(i_string))), &
922                           "- "//"Core-Core VDW, Position     [a.u.]: ", sep%fn3(i)
923                     END DO
924                  END IF
925               END SELECT
926            ELSE
927               WRITE (UNIT=output_unit, FMT="(T55,A)") "Parameters are not defined"
928            END IF
929
930            ! Additional Parameters not common to all semi-empirical methods
931            SELECT CASE (typ)
932            CASE (do_method_pdg)
933               WRITE (UNIT=output_unit, FMT="(T16,A11,T30,A,T52,F14.10,T67,F14.10)") &
934                  ADJUSTR("d_PDDG"), "- "//"Exponent [A^-1]:", dpddg/angstrom, &
935                  ADJUSTR("P_PDDG"), "- "//"Parameter  [eV]:", ppddg*evolt
936            END SELECT
937         END IF
938         CALL cp_print_key_finished_output(output_unit, logger, subsys_section, &
939                                           "PRINT%KINDS/SE_PARAMETERS")
940      END IF
941   END SUBROUTINE write_se_param
942
943END MODULE semi_empirical_types
944