LCOV - code coverage report
Current view: top level - src - pair_potential.F (source / functions) Coverage Total Hit
Test: CP2K Regtests (git:42dac4a) Lines: 97.7 % 605 591
Test Date: 2025-07-25 12:55:17 Functions: 100.0 % 7 7

            Line data    Source code
       1              : !--------------------------------------------------------------------------------------------------!
       2              : !   CP2K: A general program to perform molecular dynamics simulations                              !
       3              : !   Copyright 2000-2025 CP2K developers group <https://cp2k.org>                                   !
       4              : !                                                                                                  !
       5              : !   SPDX-License-Identifier: GPL-2.0-or-later                                                      !
       6              : !--------------------------------------------------------------------------------------------------!
       7              : 
       8              : ! **************************************************************************************************
       9              : !> \par History
      10              : !>      September 2005 - Introduced the Born-Mayer-Huggins-Fumi-Tosi  Potential (BMHTF)
      11              : !>      2006 - Major rewriting of the routines.. Linear scaling setup of splines
      12              : !>      2007 - Teodoro Laino - University of Zurich - Multiple potential
      13              : !>             Major rewriting nr.2
      14              : !> \author CJM
      15              : ! **************************************************************************************************
      16              : MODULE pair_potential
      17              : 
      18              :    USE atomic_kind_types,               ONLY: atomic_kind_type,&
      19              :                                               get_atomic_kind
      20              :    USE cp_files,                        ONLY: close_file,&
      21              :                                               open_file
      22              :    USE cp_log_handling,                 ONLY: cp_get_default_logger,&
      23              :                                               cp_logger_type,&
      24              :                                               cp_to_string
      25              :    USE fparser,                         ONLY: finalizef,&
      26              :                                               initf,&
      27              :                                               parsef
      28              :    USE kinds,                           ONLY: default_path_length,&
      29              :                                               default_string_length,&
      30              :                                               dp
      31              :    USE pair_potential_types,            ONLY: &
      32              :         ace_type, allegro_type, b4_type, bm_type, compare_pot, deepmd_type, ea_type, ft_type, &
      33              :         ftd_type, gal21_type, gal_type, gp_type, gw_type, ip_type, list_pot, lj_charmm_type, &
      34              :         lj_type, multi_type, nequip_type, nn_type, pair_potential_pp_type, &
      35              :         pair_potential_single_type, potential_single_allocation, quip_type, siepmann_type, &
      36              :         tab_type, tersoff_type, wl_type
      37              :    USE pair_potential_util,             ONLY: ener_pot,&
      38              :                                               ener_zbl,&
      39              :                                               zbl_matching_polinomial
      40              :    USE physcon,                         ONLY: bohr,&
      41              :                                               evolt,&
      42              :                                               kjmol
      43              :    USE splines_methods,                 ONLY: init_spline,&
      44              :                                               init_splinexy,&
      45              :                                               potential_s
      46              :    USE splines_types,                   ONLY: spline_data_p_type,&
      47              :                                               spline_data_type,&
      48              :                                               spline_env_create,&
      49              :                                               spline_environment_type,&
      50              :                                               spline_factor_create,&
      51              :                                               spline_factor_release,&
      52              :                                               spline_factor_type
      53              :    USE string_table,                    ONLY: str2id
      54              :    USE util,                            ONLY: sort
      55              : #include "./base/base_uses.f90"
      56              : 
      57              :    IMPLICIT NONE
      58              : 
      59              :    PRIVATE
      60              :    CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'pair_potential'
      61              :    REAL(KIND=dp), PARAMETER, PRIVATE    :: MIN_HICUT_VALUE = 1.0E-15_dp, &
      62              :                                            DEFAULT_HICUT_VALUE = 1.0E3_dp
      63              :    INTEGER, PARAMETER, PRIVATE          :: MAX_POINTS = 2000000
      64              : 
      65              :    PUBLIC :: spline_nonbond_control, &
      66              :              get_nonbond_storage, &
      67              :              init_genpot
      68              : 
      69              : CONTAINS
      70              : 
      71              : ! **************************************************************************************************
      72              : !> \brief Initialize genpot
      73              : !> \param potparm ...
      74              : !> \param ntype ...
      75              : !> \par History
      76              : !>      Teo 2007.06 - Zurich University
      77              : ! **************************************************************************************************
      78        36853 :    SUBROUTINE init_genpot(potparm, ntype)
      79              :       TYPE(pair_potential_pp_type), POINTER              :: potparm
      80              :       INTEGER, INTENT(IN)                                :: ntype
      81              : 
      82              :       CHARACTER(len=*), PARAMETER                        :: routineN = 'init_genpot'
      83              : 
      84              :       INTEGER                                            :: handle, i, j, k, ngp
      85              :       TYPE(pair_potential_single_type), POINTER          :: pot
      86              : 
      87        36853 :       CALL timeset(routineN, handle)
      88              : 
      89        36853 :       NULLIFY (pot)
      90        36853 :       ngp = 0
      91              :       ! Prescreen for general potential type
      92       896188 :       DO i = 1, ntype ! i:  first  atom type
      93     63539070 :          DO j = 1, i ! j:  second atom type
      94     62642882 :             pot => potparm%pot(i, j)%pot
      95    126145123 :             ngp = ngp + COUNT(pot%type == gp_type)
      96              :          END DO
      97              :       END DO
      98        36853 :       CALL initf(ngp)
      99        36853 :       ngp = 0
     100       896188 :       DO i = 1, ntype ! i:  first  atom type
     101     63539070 :          DO j = 1, i ! j:  second atom type
     102     62642882 :             pot => potparm%pot(i, j)%pot
     103    126145123 :             DO k = 1, SIZE(pot%type)
     104    125285788 :                IF (pot%type(k) == gp_type) THEN
     105        21972 :                   ngp = ngp + 1
     106        21972 :                   pot%set(k)%gp%myid = ngp
     107        21972 :                   CALL parsef(ngp, TRIM(pot%set(k)%gp%potential), pot%set(k)%gp%parameters)
     108              :                END IF
     109              :             END DO
     110              :          END DO
     111              :       END DO
     112        36853 :       CALL timestop(handle)
     113              : 
     114        36853 :    END SUBROUTINE init_genpot
     115              : 
     116              : ! **************************************************************************************************
     117              : !> \brief creates the splines for the potentials
     118              : !> \param spline_env ...
     119              : !> \param potparm ...
     120              : !> \param atomic_kind_set ...
     121              : !> \param eps_spline ...
     122              : !> \param max_energy ...
     123              : !> \param rlow_nb ...
     124              : !> \param emax_spline ...
     125              : !> \param npoints ...
     126              : !> \param iw ...
     127              : !> \param iw2 ...
     128              : !> \param iw3 ...
     129              : !> \param do_zbl ...
     130              : !> \param shift_cutoff ...
     131              : !> \param nonbonded_type ...
     132              : !> \par History
     133              : !>      Teo 2006.05 : Improved speed and accuracy. Linear scaling of the setup
     134              : ! **************************************************************************************************
     135         5254 :    SUBROUTINE spline_nonbond_control(spline_env, potparm, atomic_kind_set, &
     136              :                                      eps_spline, max_energy, rlow_nb, emax_spline, npoints, iw, iw2, iw3, do_zbl, &
     137              :                                      shift_cutoff, nonbonded_type)
     138              : 
     139              :       TYPE(spline_environment_type), POINTER             :: spline_env
     140              :       TYPE(pair_potential_pp_type), POINTER              :: potparm
     141              :       TYPE(atomic_kind_type), DIMENSION(:), POINTER      :: atomic_kind_set
     142              :       REAL(KIND=dp), INTENT(IN)                          :: eps_spline, max_energy, rlow_nb, &
     143              :                                                             emax_spline
     144              :       INTEGER, INTENT(IN)                                :: npoints, iw, iw2, iw3
     145              :       LOGICAL, INTENT(IN)                                :: do_zbl, shift_cutoff
     146              :       CHARACTER(LEN=*), INTENT(IN)                       :: nonbonded_type
     147              : 
     148              :       CHARACTER(len=*), PARAMETER :: routineN = 'spline_nonbond_control'
     149              : 
     150              :       INTEGER                                            :: handle, i, ip, j, k, n, ncount, &
     151              :                                                             npoints_spline, ntype
     152              :       LOGICAL                                            :: found_locut
     153              :       REAL(KIND=dp)                                      :: energy_cutoff, hicut, hicut0, locut
     154              :       TYPE(pair_potential_single_type), POINTER          :: pot
     155              : 
     156         5254 :       n = 0
     157         5254 :       ncount = 0
     158              : 
     159         5254 :       ntype = SIZE(atomic_kind_set)
     160         5254 :       CALL timeset(routineN, handle)
     161         5254 :       IF (iw3 > 0) THEN
     162              :          WRITE (iw3, "(/,T2,A,I0,A,I0,A)") &
     163         2590 :             "SPLINE_INFO| Generating ", (ntype*(ntype + 1))/2, " splines for "// &
     164         5180 :             TRIM(ADJUSTL(nonbonded_type))//" interactions "
     165              :          WRITE (iw3, "(T2,A,I0,A)") &
     166         2590 :             "             Due to ", ntype, " different atomic kinds"
     167              :       END IF
     168         5254 :       CALL init_genpot(potparm, ntype)
     169              :       ! Real computation of splines
     170         5254 :       ip = 0
     171        27664 :       DO i = 1, ntype
     172       543092 :          DO j = 1, i
     173       515428 :             pot => potparm%pot(i, j)%pot
     174       515428 :             IF (iw3 > 0 .AND. iw <= 0) THEN
     175       248590 :                IF (MOD(i*(i - 1)/2 + j, MAX(1, (ntype*(ntype + 1))/(2*10))) == 0) THEN
     176        11108 :                   WRITE (UNIT=iw3, ADVANCE="NO", FMT='(2X,A3,I0)') '...', i*(i - 1)/2 + j
     177        11108 :                   ip = ip + 1
     178        11108 :                   IF (ip >= 11) THEN
     179           96 :                      WRITE (iw3, *)
     180           96 :                      ip = 0
     181              :                   END IF
     182              :                END IF
     183              :             END IF
     184              :             ! Setup of Exclusion Types
     185       515428 :             pot%no_pp = .TRUE.
     186       515428 :             pot%no_mb = .TRUE.
     187      1030864 :             DO k = 1, SIZE(pot%type)
     188      1001147 :                SELECT CASE (pot%type(k))
     189              :                CASE (lj_type, lj_charmm_type, wl_type, gw_type, ft_type, ftd_type, ip_type, &
     190              :                      b4_type, bm_type, gp_type, ea_type, quip_type, nequip_type, allegro_type, tab_type, deepmd_type, ace_type)
     191       485711 :                   pot%no_pp = .FALSE.
     192              :                CASE (tersoff_type)
     193          118 :                   pot%no_mb = .FALSE.
     194              :                CASE (siepmann_type)
     195            5 :                   pot%no_mb = .FALSE.
     196              :                CASE (gal_type)
     197            1 :                   pot%no_mb = .FALSE.
     198              :                CASE (gal21_type)
     199            1 :                   pot%no_mb = .FALSE.
     200              :                CASE (nn_type)
     201              :                   ! Do nothing..
     202              :                CASE DEFAULT
     203              :                   ! Never reach this point
     204       515436 :                   CPABORT("")
     205              :                END SELECT
     206              :                ! Special case for EAM
     207       515428 :                SELECT CASE (pot%type(k))
     208              :                CASE (ea_type, quip_type, nequip_type, allegro_type, deepmd_type, ace_type)
     209       515436 :                   pot%no_mb = .FALSE.
     210              :                END SELECT
     211              :             END DO
     212              : 
     213              :             ! Starting SetUp of splines
     214       515428 :             IF (.NOT. pot%undef) CYCLE
     215        31565 :             ncount = ncount + 1
     216        31565 :             n = spline_env%spltab(i, j)
     217        31565 :             locut = rlow_nb
     218        31565 :             hicut0 = SQRT(pot%rcutsq)
     219        31565 :             IF (ABS(hicut0) <= MIN_HICUT_VALUE) hicut0 = DEFAULT_HICUT_VALUE
     220        31565 :             hicut = hicut0/SQRT(pot%spl_f%rcutsq_f)
     221              : 
     222        31565 :             energy_cutoff = pot%spl_f%cutoff
     223              : 
     224              :             ! Find the real locut according emax_spline
     225              :             CALL get_spline_cutoff(hicut, locut, found_locut, pot, do_zbl, &
     226        31565 :                                    energy_cutoff, emax_spline)
     227        31565 :             locut = MAX(locut*SQRT(pot%spl_f%rcutsq_f), rlow_nb)
     228              : 
     229              :             ! Real Generation of the Spline
     230        31565 :             npoints_spline = npoints
     231              :             CALL generate_spline_low(spline_env%spl_pp(n)%spl_p, npoints_spline, locut, &
     232              :                                      hicut, eps_spline, iw, iw2, i, j, n, ncount, max_energy, pot, &
     233              :                                      energy_cutoff, found_locut, do_zbl, atomic_kind_set, &
     234        31565 :                                      nonbonded_type)
     235              : 
     236        31565 :             pot%undef = .FALSE.
     237              :             ! Unique Spline working only for a pure LJ potential..
     238        31565 :             IF (SIZE(pot%type) == 1) THEN
     239        94667 :                IF (ANY(potential_single_allocation == pot%type(1))) THEN
     240              :                   ! Restoring the proper values for the generating spline pot
     241            4 :                   IF ((pot%type(1) == lj_type) .OR. (pot%type(1) == lj_charmm_type)) THEN
     242            4 :                      pot%set(1)%lj%sigma6 = pot%set(1)%lj%sigma6*pot%spl_f%rscale(1)**3
     243            4 :                      pot%set(1)%lj%sigma12 = pot%set(1)%lj%sigma6**2
     244            4 :                      pot%set(1)%lj%epsilon = pot%set(1)%lj%epsilon*pot%spl_f%fscale(1)
     245              :                   END IF
     246              :                END IF
     247              :             END IF
     248              :             ! Correct Cutoff...
     249        85540 :             IF (shift_cutoff) THEN
     250              :                pot%spl_f%cutoff = pot%spl_f%cutoff*pot%spl_f%fscale(1) - &
     251        28045 :                                   ener_pot(pot, hicut0, 0.0_dp)
     252              :             END IF
     253              :          END DO
     254              :       END DO
     255         5254 :       CALL finalizef()
     256              : 
     257         5254 :       IF (iw > 0) THEN
     258              :          WRITE (iw, '(/,T2,A,I0)') &
     259        26240 :             "SPLINE_INFO| Number of pair potential splines allocated:   ", MAXVAL(spline_env%spltab)
     260              :       END IF
     261         5254 :       IF (iw3 > 0) THEN
     262              :          WRITE (iw3, '(/,T2,A,I0,/,T2,A)') &
     263       525450 :             "SPLINE_INFO| Number of unique splines computed:            ", MAXVAL(spline_env%spltab), &
     264         5180 :             "SPLINE_INFO| Done"
     265              :       END IF
     266              : 
     267         5254 :       CALL timestop(handle)
     268              : 
     269         5254 :    END SUBROUTINE spline_nonbond_control
     270              : 
     271              : ! **************************************************************************************************
     272              : !> \brief Finds the cutoff for the generation of the spline
     273              : !>      In a two pass approach, first with low resolution, refine in a second iteration
     274              : !> \param hicut ...
     275              : !> \param locut ...
     276              : !> \param found_locut ...
     277              : !> \param pot ...
     278              : !> \param do_zbl ...
     279              : !> \param energy_cutoff ...
     280              : !> \param emax_spline ...
     281              : !> \par History
     282              : !>      Splitting in order to make some season cleaning..
     283              : !> \author Teodoro Laino [tlaino] 2007.06
     284              : ! **************************************************************************************************
     285        31565 :    SUBROUTINE get_spline_cutoff(hicut, locut, found_locut, pot, do_zbl, &
     286              :                                 energy_cutoff, emax_spline)
     287              : 
     288              :       REAL(KIND=dp), INTENT(IN)                          :: hicut
     289              :       REAL(KIND=dp), INTENT(INOUT)                       :: locut
     290              :       LOGICAL, INTENT(OUT)                               :: found_locut
     291              :       TYPE(pair_potential_single_type), OPTIONAL, &
     292              :          POINTER                                         :: pot
     293              :       LOGICAL, INTENT(IN)                                :: do_zbl
     294              :       REAL(KIND=dp), INTENT(IN)                          :: energy_cutoff, emax_spline
     295              : 
     296              :       INTEGER                                            :: ilevel, jx
     297              :       REAL(KIND=dp)                                      :: dx2, e, locut_found, x
     298              : 
     299        31565 :       dx2 = (hicut - locut)
     300        31565 :       x = hicut
     301        31565 :       locut_found = locut
     302        31565 :       found_locut = .FALSE.
     303        94695 :       DO ilevel = 1, 2
     304        63130 :          dx2 = dx2/100.0_dp
     305      5188642 :          DO jx = 1, 100
     306      5162698 :             e = ener_pot(pot, x, energy_cutoff)
     307      5162698 :             IF (do_zbl) THEN
     308         5098 :                e = e + ener_zbl(pot, x)
     309              :             END IF
     310      5162698 :             IF (ABS(e) > emax_spline) THEN
     311        37186 :                locut_found = x
     312        37186 :                found_locut = .TRUE.
     313        37186 :                EXIT
     314              :             END IF
     315      5151456 :             x = x - dx2
     316              :          END DO
     317        94695 :          x = x + dx2
     318              :       END DO
     319        31565 :       locut = locut_found
     320              : 
     321        31565 :    END SUBROUTINE get_spline_cutoff
     322              : 
     323              : ! **************************************************************************************************
     324              : !> \brief Real Generation of spline..
     325              : !> \param spl_p ...
     326              : !> \param npoints ...
     327              : !> \param locut ...
     328              : !> \param hicut ...
     329              : !> \param eps_spline ...
     330              : !> \param iw ...
     331              : !> \param iw2 ...
     332              : !> \param i ...
     333              : !> \param j ...
     334              : !> \param n ...
     335              : !> \param ncount ...
     336              : !> \param max_energy ...
     337              : !> \param pot ...
     338              : !> \param energy_cutoff ...
     339              : !> \param found_locut ...
     340              : !> \param do_zbl ...
     341              : !> \param atomic_kind_set ...
     342              : !> \param nonbonded_type ...
     343              : !> \par History
     344              : !>      Splitting in order to make some season cleaning..
     345              : !> \author Teodoro Laino [tlaino] 2007.06
     346              : ! **************************************************************************************************
     347        31565 :    SUBROUTINE generate_spline_low(spl_p, npoints, locut, hicut, eps_spline, &
     348              :                                   iw, iw2, i, j, n, ncount, max_energy, pot, energy_cutoff, &
     349              :                                   found_locut, do_zbl, atomic_kind_set, nonbonded_type)
     350              : 
     351              :       TYPE(spline_data_p_type), DIMENSION(:), POINTER    :: spl_p
     352              :       INTEGER, INTENT(INOUT)                             :: npoints
     353              :       REAL(KIND=dp), INTENT(IN)                          :: locut, hicut, eps_spline
     354              :       INTEGER, INTENT(IN)                                :: iw, iw2, i, j, n, ncount
     355              :       REAL(KIND=dp), INTENT(IN)                          :: max_energy
     356              :       TYPE(pair_potential_single_type), POINTER          :: pot
     357              :       REAL(KIND=dp), INTENT(IN), OPTIONAL                :: energy_cutoff
     358              :       LOGICAL, INTENT(IN)                                :: found_locut, do_zbl
     359              :       TYPE(atomic_kind_type), DIMENSION(:), POINTER      :: atomic_kind_set
     360              :       CHARACTER(LEN=*), INTENT(IN)                       :: nonbonded_type
     361              : 
     362              :       CHARACTER(LEN=2*default_string_length)             :: message, tmp
     363              :       CHARACTER(LEN=default_path_length)                 :: file_name
     364              :       INTEGER                                            :: ix, jx, mfac, nppa, nx, unit_number
     365              :       LOGICAL                                            :: fixed_spline_points
     366              :       REAL(KIND=dp)                                      :: df, dg, dh, diffmax, dx, dx2, e, &
     367              :                                                             e_spline, f, g, h, r, rcut, x, x2, &
     368              :                                                             xdum, xdum1, xsav
     369              :       TYPE(cp_logger_type), POINTER                      :: logger
     370              :       TYPE(spline_data_type), POINTER                    :: spline_data
     371              :       TYPE(spline_factor_type), POINTER                  :: spl_f
     372              : 
     373        31565 :       NULLIFY (logger, spl_f)
     374        63130 :       logger => cp_get_default_logger()
     375              : 
     376        31565 :       CALL spline_factor_create(spl_f)
     377        31565 :       mfac = 5
     378        31565 :       IF (npoints > 0) THEN
     379              :          fixed_spline_points = .TRUE.
     380              :       ELSE
     381        31561 :          fixed_spline_points = .FALSE.
     382        31561 :          npoints = 20
     383        31561 :          IF (.NOT. found_locut) npoints = 2
     384              :       END IF
     385        31565 :       spline_data => spl_p(1)%spline_data
     386       302539 :       DO WHILE (.TRUE.)
     387       334104 :          CALL init_splinexy(spline_data, npoints + 1)
     388       334104 :          dx2 = (1.0_dp/locut**2 - 1.0_dp/hicut**2)/REAL(npoints, KIND=dp)
     389       334104 :          x2 = 1.0_dp/hicut**2
     390       334104 :          spline_data%x1 = x2
     391    129608524 :          DO jx = 1, npoints + 1
     392              :             ! jx: loop over 1/distance**2
     393    129274420 :             x = SQRT(1.0_dp/x2)
     394    129274420 :             e = ener_pot(pot, x, energy_cutoff)
     395    129274420 :             IF (do_zbl) THEN
     396      6706340 :                e = e + ener_zbl(pot, x)
     397              :             END IF
     398    129274420 :             spline_data%y(jx) = e
     399    129608524 :             x2 = x2 + dx2
     400              :          END DO
     401       334104 :          CALL init_spline(spline_data, dx=dx2)
     402              :          ! This is the check for required accuracy on spline setup
     403       334104 :          dx2 = (hicut - locut)/REAL(mfac*npoints + 1, KIND=dp)
     404       334104 :          x2 = locut + dx2
     405       334104 :          diffmax = -1.0_dp
     406       334104 :          xsav = hicut
     407              :          ! if a fixed number of points is requested, no check on its error
     408       334104 :          IF (fixed_spline_points) EXIT
     409    644883574 :          DO jx = 1, mfac*npoints
     410    644698580 :             x = x2
     411    644698580 :             e = ener_pot(pot, x, energy_cutoff)
     412    644698580 :             IF (do_zbl) THEN
     413     33525290 :                e = e + ener_zbl(pot, x)
     414              :             END IF
     415    644698580 :             IF (ABS(e) < max_energy) THEN
     416    538831173 :                xdum1 = ABS(e - potential_s(spl_p, x*x, xdum, spl_f, logger))
     417    538831173 :                diffmax = MAX(diffmax, xdum1)
     418    538831173 :                xsav = MIN(x, xsav)
     419              :             END IF
     420    644698580 :             x2 = x2 + dx2
     421    644883574 :             IF (x2 > hicut) EXIT
     422              :          END DO
     423       334100 :          IF (npoints > MAX_POINTS) THEN
     424            0 :             WRITE (message, '(A,I8,A,G12.6,A)') "SPLINE_INFO| Number of points: ", npoints, &
     425            0 :                " obtained accuracy ", diffmax, ". MM SPLINE: no convergence on required"// &
     426            0 :                " accuracy (adjust EPS_SPLINE and rerun)"
     427            0 :             CALL cp_abort(__LOCATION__, TRIM(message))
     428              :          END IF
     429              :          ! accuracy is poor or we have found no points below max_energy, refine mesh
     430       334104 :          IF (diffmax > eps_spline .OR. diffmax < 0.0_dp) THEN
     431       302539 :             npoints = CEILING(1.2_dp*REAL(npoints, KIND=dp))
     432              :          ELSE
     433              :             EXIT
     434              :          END IF
     435              :       END DO
     436              :       ! Print spline info to STDOUT if requested
     437        31565 :       IF (iw > 0) THEN
     438              :          WRITE (UNIT=iw, &
     439              :                 FMT="(/,A,I0,/,A,I0,/,A,I0,1X,I0,/,A,/,A,I0,2(/,A,ES13.6),2(/,A,2ES13.6))") &
     440         4800 :             " SPLINE_INFO| Spline number:                                ", ncount, &
     441         4800 :             " SPLINE_INFO| Unique spline number:                         ", n, &
     442         4800 :             " SPLINE_INFO| Atomic kind numbers:                          ", i, j, &
     443              :             " SPLINE_INFO| Atomic kind names:                            "//TRIM(ADJUSTL(atomic_kind_set(i)%name))//" "// &
     444         4800 :             TRIM(ADJUSTL(atomic_kind_set(j)%name)), &
     445         4800 :             " SPLINE_INFO| Number of spline points:                      ", npoints, &
     446         4800 :             " SPLINE_INFO| Requested accuracy [Hartree]:                ", eps_spline, &
     447         4800 :             " SPLINE_INFO| Achieved accuracy [Hartree]:                 ", diffmax, &
     448         4800 :             " SPLINE_INFO| Spline range [bohr]:                         ", locut, hicut, &
     449         9600 :             " SPLINE_INFO| Spline range used to achieve accuracy [bohr]:", xsav, hicut
     450         4800 :          dx2 = (hicut - locut)/REAL(npoints + 1, KIND=dp)
     451         4800 :          x = locut + dx2
     452              :          WRITE (UNIT=iw, FMT='(A,ES17.9)') &
     453         4800 :             " SPLINE_INFO| Spline value at RMIN [Hartree]:             ", potential_s(spl_p, x*x, xdum, spl_f, logger), &
     454         4800 :             " SPLINE_INFO| Spline value at RMAX [Hartree]:             ", potential_s(spl_p, hicut*hicut, xdum, spl_f, logger), &
     455         9600 :             " SPLINE_INFO| Non-bonded energy cutoff [Hartree]:         ", energy_cutoff
     456              :       END IF
     457              :       ! Print spline data on file if requested
     458        31565 :       IF (iw2 > 0) THEN
     459              :          ! Set increment to 200 points per Angstrom
     460           64 :          nppa = 200
     461           64 :          dx = bohr/REAL(nppa, KIND=dp)
     462           64 :          nx = NINT(hicut/dx)
     463           64 :          file_name = ""
     464           64 :          tmp = ADJUSTL(cp_to_string(n))
     465              :          WRITE (UNIT=file_name, FMT="(A,I0,A)") &
     466              :             TRIM(ADJUSTL(nonbonded_type))//"_SPLINE_"//TRIM(tmp)//"_"// &
     467              :             TRIM(ADJUSTL(atomic_kind_set(i)%name))//"_"// &
     468           64 :             TRIM(ADJUSTL(atomic_kind_set(j)%name))
     469              :          CALL open_file(file_name=file_name, &
     470              :                         file_status="UNKNOWN", &
     471              :                         file_form="FORMATTED", &
     472              :                         file_action="WRITE", &
     473           64 :                         unit_number=unit_number)
     474              :          WRITE (UNIT=unit_number, &
     475              :                 FMT="(2(A,I0,/),A,I0,1X,I0,/,A,/,A,I0,2(/,A,ES13.6),2(/,A,2ES13.6),/,A,ES13.6,/,A,I0,A,/,A)") &
     476           64 :             "# Spline number:                                    ", ncount, &
     477           64 :             "# Unique spline number:                             ", n, &
     478           64 :             "# Atomic kind numbers:                              ", i, j, &
     479              :             "# Atomic kind names:                                "//TRIM(ADJUSTL(atomic_kind_set(i)%name))//" "// &
     480           64 :             TRIM(ADJUSTL(atomic_kind_set(j)%name)), &
     481           64 :             "# Number of spline points:                          ", npoints, &
     482           64 :             "# Requested accuracy [eV]:                         ", eps_spline*evolt, &
     483           64 :             "# Achieved accuracy [eV]:                          ", diffmax*evolt, &
     484           64 :             "# Spline range [Angstrom]:                         ", locut/bohr, hicut/bohr, &
     485           64 :             "# Spline range used to achieve accuracy [Angstrom]:", xsav/bohr, hicut/bohr, &
     486           64 :             "# Non-bonded energy cutoff [eV]:                   ", energy_cutoff*evolt, &
     487           64 :             "# Test spline using ", nppa, " points per Angstrom:", &
     488              :             "#     Abscissa [Angstrom]              Energy [eV]      Splined energy [eV] Derivative [eV/Angstrom]"// &
     489          128 :             "      |Energy error| [eV]"
     490           64 :          x = 0.0_dp
     491       128296 :          DO jx = 0, nx
     492       128232 :             IF (x > hicut) x = hicut
     493       128232 :             IF (x > locut) THEN
     494       106526 :                e = ener_pot(pot, x, energy_cutoff)
     495       106526 :                IF (do_zbl) e = e + ener_zbl(pot, x)
     496       106526 :                e_spline = potential_s(spl_p, x*x, xdum, spl_f, logger)
     497              :                WRITE (UNIT=unit_number, FMT="(5ES25.12)") &
     498       106526 :                   x/bohr, e*evolt, e_spline*evolt, -bohr*x*xdum*evolt, ABS((e - e_spline)*evolt)
     499              :             END IF
     500       128296 :             x = x + dx
     501              :          END DO
     502           64 :          CALL close_file(unit_number=unit_number)
     503              :          !MK Write table.xvf for GROMACS 4.5.5
     504              :          WRITE (UNIT=file_name, FMT="(A,I0,A)") &
     505              :             "table_"// &
     506              :             TRIM(ADJUSTL(atomic_kind_set(i)%name))//"_"// &
     507           64 :             TRIM(ADJUSTL(atomic_kind_set(j)%name))//".xvg"
     508              :          CALL open_file(file_name=file_name, &
     509              :                         file_status="UNKNOWN", &
     510              :                         file_form="FORMATTED", &
     511              :                         file_action="WRITE", &
     512           64 :                         unit_number=unit_number)
     513              :          ! Recommended increment for dp is 0.0005 nm = 0.005 Angstrom
     514              :          ! which are 200 points/Angstrom
     515           64 :          rcut = 0.1_dp*hicut/bohr
     516           64 :          x = 0.0_dp
     517       128296 :          DO jx = 0, nx
     518       128232 :             IF (x > hicut) x = hicut
     519       128232 :             r = 0.1_dp*x/bohr ! Convert bohr to nm
     520       128232 :             IF (x <= locut) THEN
     521              :                WRITE (UNIT=unit_number, FMT="(7ES25.12)") &
     522       151942 :                   r, (0.0_dp, ix=1, 6)
     523              :             ELSE
     524       106526 :                e_spline = potential_s(spl_p, x*x, xdum, spl_f, logger)
     525       106526 :                f = 1.0_dp/r
     526       106526 :                df = -1.0_dp/r**2
     527       106526 :                g = -1.0_dp/r**6 + 1.0_dp/rcut**6
     528       106526 :                dg = 6.0_dp/r**7
     529       106526 :                h = e_spline*kjmol
     530       106526 :                dh = -10.0_dp*bohr*x*xdum*kjmol
     531              :                WRITE (UNIT=unit_number, FMT="(7ES25.12)") &
     532       106526 :                   r, f, -df, & ! r, f(r), -f'(r) => probably not used
     533       106526 :                   g, -dg, & !    g(r), -g'(r) => not used, if C = 0
     534       213052 :                   h, -dh !    h(r), -h'(r) => used, if A = 1
     535              :             END IF
     536       128296 :             x = x + dx
     537              :          END DO
     538           64 :          CALL close_file(unit_number=unit_number)
     539              :       END IF
     540              : 
     541        31565 :       CALL spline_factor_release(spl_f)
     542              : 
     543        31565 :    END SUBROUTINE generate_spline_low
     544              : 
     545              : ! **************************************************************************************************
     546              : !> \brief Prescreening of the effective bonds evaluations. linear scaling algorithm
     547              : !> \param spline_env ...
     548              : !> \param potparm ...
     549              : !> \param atomic_kind_set ...
     550              : !> \param do_zbl ...
     551              : !> \param shift_cutoff ...
     552              : !> \author Teodoro Laino [tlaino] 2006.05
     553              : ! **************************************************************************************************
     554         5254 :    SUBROUTINE get_nonbond_storage(spline_env, potparm, atomic_kind_set, do_zbl, &
     555              :                                   shift_cutoff)
     556              : 
     557              :       TYPE(spline_environment_type), POINTER             :: spline_env
     558              :       TYPE(pair_potential_pp_type), POINTER              :: potparm
     559              :       TYPE(atomic_kind_type), DIMENSION(:), POINTER      :: atomic_kind_set
     560              :       LOGICAL, INTENT(IN)                                :: do_zbl, shift_cutoff
     561              : 
     562              :       CHARACTER(len=*), PARAMETER :: routineN = 'get_nonbond_storage'
     563              : 
     564              :       INTEGER                                            :: handle, i, idim, iend, istart, j, k, &
     565              :                                                             locij, n, ndim, nk, ntype, nunique, &
     566              :                                                             nvar, pot_target, tmpij(2), tmpij0(2)
     567         5254 :       INTEGER, ALLOCATABLE, DIMENSION(:)                 :: Iwork1, Iwork2, my_index
     568         5254 :       INTEGER, ALLOCATABLE, DIMENSION(:, :)              :: tmp_index
     569              :       LOGICAL                                            :: at_least_one, check
     570         5254 :       REAL(KIND=dp), ALLOCATABLE, DIMENSION(:)           :: Cwork, Rwork, wtmp
     571         5254 :       REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :)        :: pot_par
     572              : 
     573         5254 :       CALL timeset(routineN, handle)
     574              : 
     575         5254 :       ntype = SIZE(atomic_kind_set)
     576        27664 :       DO i = 1, ntype
     577       543092 :          DO j = 1, i
     578       537838 :             potparm%pot(i, j)%pot%undef = .FALSE.
     579              :          END DO
     580              :       END DO
     581        21016 :       ALLOCATE (tmp_index(ntype, ntype))
     582              :       !
     583         5254 :       nunique = 0
     584      1036110 :       tmp_index = HUGE(0)
     585       120842 :       DO pot_target = MINVAL(list_pot), MAXVAL(list_pot)
     586       115588 :          ndim = 0
     587       608608 :          DO i = 1, ntype
     588     11948024 :             DO j = 1, i
     589     11339416 :                IF (SIZE(potparm%pot(i, j)%pot%type) /= 1) CYCLE
     590     11832260 :                IF (potparm%pot(i, j)%pot%type(1) == pot_target) THEN
     591       515420 :                   tmp_index(i, j) = 1
     592       515420 :                   tmp_index(j, i) = 1
     593       515420 :                   ndim = ndim + 1
     594              :                END IF
     595              :             END DO
     596              :          END DO
     597       115588 :          IF (ndim == 0) CYCLE ! No potential of this kind found
     598         6159 :          nvar = 0
     599              :          SELECT CASE (pot_target)
     600              :          CASE (lj_type, lj_charmm_type)
     601              :             nvar = 3 + nvar
     602              :          CASE (wl_type)
     603            0 :             nvar = 3 + nvar
     604              :          CASE (gw_type)
     605            0 :             nvar = 5 + nvar
     606              :          CASE (ea_type)
     607           12 :             nvar = 4 + nvar
     608              :          CASE (quip_type)
     609            0 :             nvar = 1 + nvar
     610              :          CASE (nequip_type)
     611            4 :             nvar = 1 + nvar
     612              :          CASE (allegro_type)
     613            4 :             nvar = 1 + nvar
     614              :          CASE (ace_type)
     615            6 :             nvar = 2 + nvar
     616              :          CASE (deepmd_type)
     617            2 :             nvar = 2 + nvar
     618              :          CASE (ft_type)
     619            4 :             nvar = 4 + nvar
     620              :          CASE (ftd_type)
     621           18 :             nvar = 6 + nvar
     622              :          CASE (ip_type)
     623          260 :             nvar = 3 + nvar
     624              :          CASE (b4_type)
     625          260 :             nvar = 6 + nvar
     626              :          CASE (bm_type)
     627            6 :             nvar = 9 + nvar
     628              :          CASE (gp_type)
     629          568 :             nvar = 2 + nvar
     630              :          CASE (tersoff_type)
     631           38 :             nvar = 13 + nvar
     632              :          CASE (siepmann_type)
     633            5 :             nvar = 5 + nvar
     634              :          CASE (gal_type)
     635            1 :             nvar = 12 + nvar
     636              :          CASE (gal21_type)
     637            1 :             nvar = 30 + nvar
     638              :          CASE (nn_type)
     639         2063 :             nvar = nvar
     640              :          CASE (tab_type)
     641            8 :             nvar = 4 + nvar
     642              :          CASE DEFAULT
     643         6159 :             CPABORT("")
     644              :          END SELECT
     645              :          ! Setup a table of the indexes..
     646        18477 :          ALLOCATE (my_index(ndim))
     647         6159 :          n = 0
     648         6159 :          nk = 0
     649        35760 :          DO i = 1, ntype
     650       971850 :             DO j = 1, i
     651       936090 :                n = n + 1
     652       936090 :                IF (SIZE(potparm%pot(i, j)%pot%type) /= 1) CYCLE
     653       965687 :                IF (potparm%pot(i, j)%pot%type(1) == pot_target) THEN
     654       515420 :                   nk = nk + 1
     655       515420 :                   my_index(nk) = n
     656              :                END IF
     657              :             END DO
     658              :          END DO
     659         6159 :          IF (nvar /= 0) THEN
     660        16384 :             ALLOCATE (pot_par(ndim, nvar))
     661         4096 :             n = 0
     662         4096 :             nk = 0
     663        23656 :             DO i = 1, ntype
     664       532273 :                DO j = 1, i
     665       508617 :                   n = n + 1
     666       508617 :                   IF (SIZE(potparm%pot(i, j)%pot%type) /= 1) CYCLE
     667       528173 :                   IF (potparm%pot(i, j)%pot%type(1) == pot_target) THEN
     668       485820 :                      nk = nk + 1
     669       485820 :                      my_index(nk) = n
     670       480956 :                      SELECT CASE (pot_target)
     671              :                      CASE (lj_type, lj_charmm_type)
     672       480956 :                         pot_par(nk, 1) = potparm%pot(i, j)%pot%set(1)%lj%epsilon
     673       480956 :                         pot_par(nk, 2) = potparm%pot(i, j)%pot%set(1)%lj%sigma6
     674       480956 :                         pot_par(nk, 3) = potparm%pot(i, j)%pot%set(1)%lj%sigma12
     675              :                      CASE (gp_type)
     676         3208 :                         pot_par(nk, 1) = str2id(potparm%pot(i, j)%pot%set(1)%gp%potential)
     677         3208 :                         pot_par(nk, 2) = str2id(potparm%pot(i, j)%pot%set(1)%gp%variables)
     678              :                      CASE (wl_type)
     679         1049 :                         pot_par(nk, 1) = potparm%pot(i, j)%pot%set(1)%willis%a
     680         1049 :                         pot_par(nk, 2) = potparm%pot(i, j)%pot%set(1)%willis%b
     681         1049 :                         pot_par(nk, 3) = potparm%pot(i, j)%pot%set(1)%willis%c
     682              :                      CASE (gw_type)
     683            0 :                         pot_par(nk, 1) = potparm%pot(i, j)%pot%set(1)%goodwin%vr0
     684            0 :                         pot_par(nk, 2) = potparm%pot(i, j)%pot%set(1)%goodwin%m
     685            0 :                         pot_par(nk, 3) = potparm%pot(i, j)%pot%set(1)%goodwin%mc
     686            0 :                         pot_par(nk, 4) = potparm%pot(i, j)%pot%set(1)%goodwin%d
     687            0 :                         pot_par(nk, 5) = potparm%pot(i, j)%pot%set(1)%goodwin%dc
     688              :                      CASE (ea_type)
     689           20 :                         pot_par(nk, 1) = potparm%pot(i, j)%pot%set(1)%eam%drar
     690           20 :                         pot_par(nk, 2) = potparm%pot(i, j)%pot%set(1)%eam%drhoar
     691           20 :                         pot_par(nk, 3) = potparm%pot(i, j)%pot%set(1)%eam%acutal
     692           20 :                         pot_par(nk, 4) = potparm%pot(i, j)%pot%set(1)%eam%npoints
     693              :                      CASE (quip_type)
     694              :                         pot_par(nk, 1) = str2id( &
     695              :                                          TRIM(potparm%pot(i, j)%pot%set(1)%quip%quip_file_name)// &
     696              :                                          TRIM(potparm%pot(i, j)%pot%set(1)%quip%init_args)// &
     697            0 :                                          TRIM(potparm%pot(i, j)%pot%set(1)%quip%calc_args))
     698              :                      CASE (nequip_type)
     699              :                         pot_par(nk, 1) = str2id( &
     700           12 :                                          TRIM(potparm%pot(i, j)%pot%set(1)%nequip%nequip_file_name))
     701              :                      CASE (allegro_type)
     702              :                         pot_par(nk, 1) = str2id( &
     703            8 :                                          TRIM(potparm%pot(i, j)%pot%set(1)%allegro%allegro_file_name))
     704              :                      CASE (ace_type)
     705              :                         pot_par(nk, 1) = str2id( &
     706           18 :                                          TRIM(potparm%pot(i, j)%pot%set(1)%ace%ace_file_name))
     707           18 :                         pot_par(nk, 2) = potparm%pot(i, j)%pot%set(1)%ace%atom_ace_type
     708              :                      CASE (deepmd_type)
     709              :                         pot_par(nk, 1) = str2id( &
     710            6 :                                          TRIM(potparm%pot(i, j)%pot%set(1)%deepmd%deepmd_file_name))
     711            6 :                         pot_par(nk, 2) = potparm%pot(i, j)%pot%set(1)%deepmd%atom_deepmd_type
     712              :                      CASE (ft_type)
     713           12 :                         pot_par(nk, 1) = potparm%pot(i, j)%pot%set(1)%ft%A
     714           12 :                         pot_par(nk, 2) = potparm%pot(i, j)%pot%set(1)%ft%B
     715           12 :                         pot_par(nk, 3) = potparm%pot(i, j)%pot%set(1)%ft%C
     716           12 :                         pot_par(nk, 4) = potparm%pot(i, j)%pot%set(1)%ft%D
     717              :                      CASE (ftd_type)
     718           66 :                         pot_par(nk, 1) = potparm%pot(i, j)%pot%set(1)%ftd%A
     719           66 :                         pot_par(nk, 2) = potparm%pot(i, j)%pot%set(1)%ftd%B
     720           66 :                         pot_par(nk, 3) = potparm%pot(i, j)%pot%set(1)%ftd%C
     721           66 :                         pot_par(nk, 4) = potparm%pot(i, j)%pot%set(1)%ftd%D
     722           66 :                         pot_par(nk, 5) = potparm%pot(i, j)%pot%set(1)%ftd%BD(1)
     723           66 :                         pot_par(nk, 6) = potparm%pot(i, j)%pot%set(1)%ftd%BD(2)
     724              :                      CASE (ip_type)
     725           48 :                         pot_par(nk, 1) = potparm%pot(i, j)%pot%set(1)%ipbv%rcore
     726           48 :                         pot_par(nk, 2) = potparm%pot(i, j)%pot%set(1)%ipbv%m
     727           48 :                         pot_par(nk, 3) = potparm%pot(i, j)%pot%set(1)%ipbv%b
     728              :                      CASE (b4_type)
     729          260 :                         pot_par(nk, 1) = potparm%pot(i, j)%pot%set(1)%buck4r%a
     730          260 :                         pot_par(nk, 2) = potparm%pot(i, j)%pot%set(1)%buck4r%b
     731          260 :                         pot_par(nk, 3) = potparm%pot(i, j)%pot%set(1)%buck4r%c
     732          260 :                         pot_par(nk, 4) = potparm%pot(i, j)%pot%set(1)%buck4r%r1
     733          260 :                         pot_par(nk, 5) = potparm%pot(i, j)%pot%set(1)%buck4r%r2
     734          260 :                         pot_par(nk, 6) = potparm%pot(i, j)%pot%set(1)%buck4r%r3
     735              :                      CASE (bm_type)
     736           10 :                         pot_par(nk, 1) = potparm%pot(i, j)%pot%set(1)%buckmo%f0
     737           10 :                         pot_par(nk, 2) = potparm%pot(i, j)%pot%set(1)%buckmo%a1
     738           10 :                         pot_par(nk, 3) = potparm%pot(i, j)%pot%set(1)%buckmo%a2
     739           10 :                         pot_par(nk, 4) = potparm%pot(i, j)%pot%set(1)%buckmo%b1
     740           10 :                         pot_par(nk, 5) = potparm%pot(i, j)%pot%set(1)%buckmo%b2
     741           10 :                         pot_par(nk, 6) = potparm%pot(i, j)%pot%set(1)%buckmo%c
     742           10 :                         pot_par(nk, 7) = potparm%pot(i, j)%pot%set(1)%buckmo%d
     743           10 :                         pot_par(nk, 8) = potparm%pot(i, j)%pot%set(1)%buckmo%r0
     744           10 :                         pot_par(nk, 9) = potparm%pot(i, j)%pot%set(1)%buckmo%beta
     745              :                      CASE (tersoff_type)
     746          116 :                         pot_par(nk, 1) = potparm%pot(i, j)%pot%set(1)%tersoff%A
     747          116 :                         pot_par(nk, 2) = potparm%pot(i, j)%pot%set(1)%tersoff%B
     748          116 :                         pot_par(nk, 3) = potparm%pot(i, j)%pot%set(1)%tersoff%lambda1
     749          116 :                         pot_par(nk, 4) = potparm%pot(i, j)%pot%set(1)%tersoff%lambda2
     750          116 :                         pot_par(nk, 5) = potparm%pot(i, j)%pot%set(1)%tersoff%alpha
     751          116 :                         pot_par(nk, 6) = potparm%pot(i, j)%pot%set(1)%tersoff%beta
     752          116 :                         pot_par(nk, 7) = potparm%pot(i, j)%pot%set(1)%tersoff%n
     753          116 :                         pot_par(nk, 8) = potparm%pot(i, j)%pot%set(1)%tersoff%c
     754          116 :                         pot_par(nk, 9) = potparm%pot(i, j)%pot%set(1)%tersoff%d
     755          116 :                         pot_par(nk, 10) = potparm%pot(i, j)%pot%set(1)%tersoff%h
     756          116 :                         pot_par(nk, 11) = potparm%pot(i, j)%pot%set(1)%tersoff%lambda3
     757          116 :                         pot_par(nk, 12) = potparm%pot(i, j)%pot%set(1)%tersoff%bigR
     758          116 :                         pot_par(nk, 13) = potparm%pot(i, j)%pot%set(1)%tersoff%bigD
     759              :                      CASE (siepmann_type)
     760            5 :                         pot_par(nk, 1) = potparm%pot(i, j)%pot%set(1)%siepmann%B
     761            5 :                         pot_par(nk, 2) = potparm%pot(i, j)%pot%set(1)%siepmann%D
     762            5 :                         pot_par(nk, 3) = potparm%pot(i, j)%pot%set(1)%siepmann%E
     763            5 :                         pot_par(nk, 4) = potparm%pot(i, j)%pot%set(1)%siepmann%F
     764            5 :                         pot_par(nk, 5) = potparm%pot(i, j)%pot%set(1)%siepmann%beta
     765              :                      CASE (gal_type)
     766            1 :                         pot_par(nk, 1) = potparm%pot(i, j)%pot%set(1)%gal%epsilon
     767            1 :                         pot_par(nk, 2) = potparm%pot(i, j)%pot%set(1)%gal%bxy
     768            1 :                         pot_par(nk, 3) = potparm%pot(i, j)%pot%set(1)%gal%bz
     769            1 :                         pot_par(nk, 4) = potparm%pot(i, j)%pot%set(1)%gal%r1
     770            1 :                         pot_par(nk, 5) = potparm%pot(i, j)%pot%set(1)%gal%r2
     771            1 :                         pot_par(nk, 6) = potparm%pot(i, j)%pot%set(1)%gal%a1
     772            1 :                         pot_par(nk, 7) = potparm%pot(i, j)%pot%set(1)%gal%a2
     773            1 :                         pot_par(nk, 8) = potparm%pot(i, j)%pot%set(1)%gal%a3
     774            1 :                         pot_par(nk, 9) = potparm%pot(i, j)%pot%set(1)%gal%a4
     775            1 :                         pot_par(nk, 10) = potparm%pot(i, j)%pot%set(1)%gal%a
     776            1 :                         pot_par(nk, 11) = potparm%pot(i, j)%pot%set(1)%gal%b
     777            1 :                         pot_par(nk, 12) = potparm%pot(i, j)%pot%set(1)%gal%c
     778              :                      CASE (gal21_type)
     779            1 :                         pot_par(nk, 1) = potparm%pot(i, j)%pot%set(1)%gal21%epsilon1
     780            1 :                         pot_par(nk, 2) = potparm%pot(i, j)%pot%set(1)%gal21%epsilon2
     781            1 :                         pot_par(nk, 3) = potparm%pot(i, j)%pot%set(1)%gal21%epsilon3
     782            1 :                         pot_par(nk, 4) = potparm%pot(i, j)%pot%set(1)%gal21%bxy1
     783            1 :                         pot_par(nk, 5) = potparm%pot(i, j)%pot%set(1)%gal21%bxy2
     784            1 :                         pot_par(nk, 6) = potparm%pot(i, j)%pot%set(1)%gal21%bz1
     785            1 :                         pot_par(nk, 7) = potparm%pot(i, j)%pot%set(1)%gal21%bz2
     786            1 :                         pot_par(nk, 8) = potparm%pot(i, j)%pot%set(1)%gal21%r1
     787            1 :                         pot_par(nk, 9) = potparm%pot(i, j)%pot%set(1)%gal21%r2
     788            1 :                         pot_par(nk, 10) = potparm%pot(i, j)%pot%set(1)%gal21%a11
     789            1 :                         pot_par(nk, 11) = potparm%pot(i, j)%pot%set(1)%gal21%a12
     790            1 :                         pot_par(nk, 12) = potparm%pot(i, j)%pot%set(1)%gal21%a13
     791            1 :                         pot_par(nk, 13) = potparm%pot(i, j)%pot%set(1)%gal21%a21
     792            1 :                         pot_par(nk, 14) = potparm%pot(i, j)%pot%set(1)%gal21%a22
     793            1 :                         pot_par(nk, 15) = potparm%pot(i, j)%pot%set(1)%gal21%a23
     794            1 :                         pot_par(nk, 16) = potparm%pot(i, j)%pot%set(1)%gal21%a31
     795            1 :                         pot_par(nk, 17) = potparm%pot(i, j)%pot%set(1)%gal21%a32
     796            1 :                         pot_par(nk, 18) = potparm%pot(i, j)%pot%set(1)%gal21%a33
     797            1 :                         pot_par(nk, 19) = potparm%pot(i, j)%pot%set(1)%gal21%a41
     798            1 :                         pot_par(nk, 20) = potparm%pot(i, j)%pot%set(1)%gal21%a42
     799            1 :                         pot_par(nk, 21) = potparm%pot(i, j)%pot%set(1)%gal21%a43
     800            1 :                         pot_par(nk, 22) = potparm%pot(i, j)%pot%set(1)%gal21%AO1
     801            1 :                         pot_par(nk, 23) = potparm%pot(i, j)%pot%set(1)%gal21%AO2
     802            1 :                         pot_par(nk, 24) = potparm%pot(i, j)%pot%set(1)%gal21%BO1
     803            1 :                         pot_par(nk, 25) = potparm%pot(i, j)%pot%set(1)%gal21%BO2
     804            1 :                         pot_par(nk, 26) = potparm%pot(i, j)%pot%set(1)%gal21%c
     805            1 :                         pot_par(nk, 27) = potparm%pot(i, j)%pot%set(1)%gal21%AH1
     806            1 :                         pot_par(nk, 28) = potparm%pot(i, j)%pot%set(1)%gal21%AH2
     807            1 :                         pot_par(nk, 29) = potparm%pot(i, j)%pot%set(1)%gal21%BH1
     808            1 :                         pot_par(nk, 30) = potparm%pot(i, j)%pot%set(1)%gal21%BH2
     809              :                      CASE (tab_type)
     810           24 :                         pot_par(nk, 1) = potparm%pot(i, j)%pot%set(1)%tab%dr
     811           24 :                         pot_par(nk, 2) = potparm%pot(i, j)%pot%set(1)%tab%rcut
     812           24 :                         pot_par(nk, 3) = potparm%pot(i, j)%pot%set(1)%tab%npoints
     813           24 :                         pot_par(nk, 4) = potparm%pot(i, j)%pot%set(1)%tab%index
     814              :                      CASE (nn_type)
     815              :                         ! no checks
     816              :                      CASE DEFAULT
     817       485820 :                         CPABORT("")
     818              :                      END SELECT
     819      1448076 :                      IF (ANY(potential_single_allocation == pot_target)) THEN
     820        37536 :                         pot_par(nk, :) = REAL(pot_target, KIND=dp)
     821              :                      END IF
     822              :                   END IF
     823              :                END DO
     824              :             END DO
     825              :             ! Main Sorting Loop
     826        12288 :             ALLOCATE (Rwork(ndim))
     827         8192 :             ALLOCATE (Iwork1(ndim))
     828         8192 :             ALLOCATE (Iwork2(ndim))
     829         8192 :             ALLOCATE (wtmp(nvar))
     830         4096 :             CALL sort(pot_par(:, 1), ndim, Iwork1)
     831              :             ! Sort all the other components of the potential
     832        13016 :             DO k = 2, nvar
     833       979588 :                Rwork(:) = pot_par(:, k)
     834       983684 :                DO i = 1, ndim
     835       979588 :                   pot_par(i, k) = Rwork(Iwork1(i))
     836              :                END DO
     837              :             END DO
     838       489916 :             Iwork2(:) = my_index
     839       489916 :             DO i = 1, ndim
     840       489916 :                my_index(i) = Iwork2(Iwork1(i))
     841              :             END DO
     842              :             ! Iterative sorting
     843         7695 :             DO k = 2, nvar
     844        13149 :                wtmp(1:k - 1) = pot_par(1, 1:k - 1)
     845              :                istart = 1
     846              :                at_least_one = .FALSE.
     847       969920 :                DO j = 1, ndim
     848       964215 :                   Rwork(j) = pot_par(j, k)
     849      2362854 :                   IF (ALL(pot_par(j, 1:k - 1) == wtmp(1:k - 1))) CYCLE
     850        35506 :                   iend = j - 1
     851        90797 :                   wtmp(1:k - 1) = pot_par(j, 1:k - 1)
     852              :                   ! If the ordered array has no two same consecutive elements
     853              :                   ! does not make any sense to proceed ordering the others
     854              :                   ! related parameters..
     855        35506 :                   idim = iend - istart + 1
     856        35506 :                   CALL sort(Rwork(istart:iend), idim, Iwork1(istart:iend))
     857       967420 :                   Iwork1(istart:iend) = Iwork1(istart:iend) - 1 + istart
     858        35506 :                   IF (idim /= 1) at_least_one = .TRUE.
     859       934414 :                   istart = j
     860              :                END DO
     861         5705 :                iend = ndim
     862         5705 :                idim = iend - istart + 1
     863         5705 :                CALL sort(Rwork(istart:iend), idim, Iwork1(istart:iend))
     864        38006 :                Iwork1(istart:iend) = Iwork1(istart:iend) - 1 + istart
     865         5705 :                IF (idim /= 1) at_least_one = .TRUE.
     866       969920 :                pot_par(:, k) = Rwork
     867         5705 :                IF (.NOT. at_least_one) EXIT
     868              :                ! Sort other components
     869         5358 :                DO j = k + 1, nvar
     870       484450 :                   Rwork(:) = pot_par(:, j)
     871       488049 :                   DO i = 1, ndim
     872       484450 :                      pot_par(i, j) = Rwork(Iwork1(i))
     873              :                   END DO
     874              :                END DO
     875       962282 :                Iwork2(:) = my_index
     876       966378 :                DO i = 1, ndim
     877       962282 :                   my_index(i) = Iwork2(Iwork1(i))
     878              :                END DO
     879              :             END DO
     880         4096 :             DEALLOCATE (wtmp)
     881         4096 :             DEALLOCATE (Iwork1)
     882         4096 :             DEALLOCATE (Iwork2)
     883         4096 :             DEALLOCATE (Rwork)
     884              :             !
     885              :             ! Let's determine the number of unique potentials and tag them
     886              :             !
     887         8192 :             ALLOCATE (Cwork(nvar))
     888        17112 :             Cwork(:) = pot_par(1, :)
     889         4096 :             locij = my_index(1)
     890         4096 :             CALL get_indexes(locij, ntype, tmpij0)
     891         4096 :             istart = 1
     892       489916 :             DO j = 1, ndim
     893              :                ! Special cases for EAM and IPBV
     894       485820 :                locij = my_index(j)
     895       485820 :                CALL get_indexes(locij, ntype, tmpij)
     896           68 :                SELECT CASE (pot_target)
     897              :                   !NB should do something about QUIP here?
     898              :                CASE (ea_type, ip_type)
     899              :                   ! check the array components
     900              :                   CALL compare_pot(potparm%pot(tmpij(1), tmpij(2))%pot, &
     901              :                                    potparm%pot(tmpij0(1), tmpij0(2))%pot, &
     902         3276 :                                    check)
     903              :                CASE (gp_type)
     904         3208 :                   check = .TRUE.
     905         3208 :                   IF (ASSOCIATED(potparm%pot(tmpij(1), tmpij(2))%pot%set(1)%gp%parameters) .AND. &
     906              :                       ASSOCIATED(potparm%pot(tmpij0(1), tmpij0(2))%pot%set(1)%gp%parameters)) THEN
     907         3208 :                      IF (SIZE(potparm%pot(tmpij(1), tmpij(2))%pot%set(1)%gp%parameters) == &
     908              :                          SIZE(potparm%pot(tmpij0(1), tmpij0(2))%pot%set(1)%gp%parameters)) THEN
     909        12640 :                         IF (ANY(potparm%pot(tmpij(1), tmpij(2))%pot%set(1)%gp%parameters /= &
     910            0 :                                 potparm%pot(tmpij0(1), tmpij0(2))%pot%set(1)%gp%parameters)) check = .FALSE.
     911              :                      END IF
     912              :                   END IF
     913         3208 :                   IF (ASSOCIATED(potparm%pot(tmpij(1), tmpij(2))%pot%set(1)%gp%values) .AND. &
     914       482544 :                       ASSOCIATED(potparm%pot(tmpij0(1), tmpij0(2))%pot%set(1)%gp%values)) THEN
     915         3208 :                      IF (SIZE(potparm%pot(tmpij(1), tmpij(2))%pot%set(1)%gp%values) == &
     916              :                          SIZE(potparm%pot(tmpij0(1), tmpij0(2))%pot%set(1)%gp%values)) THEN
     917         7502 :                         IF (ANY(potparm%pot(tmpij(1), tmpij(2))%pot%set(1)%gp%values /= &
     918         2569 :                                 potparm%pot(tmpij0(1), tmpij0(2))%pot%set(1)%gp%values)) check = .FALSE.
     919              :                      END IF
     920              :                   END IF
     921              :                CASE default
     922       485820 :                   check = .TRUE.
     923              :                END SELECT
     924      1880737 :                IF (ALL(Cwork == pot_par(j, :)) .AND. check) CYCLE
     925        99223 :                Cwork(:) = pot_par(j, :)
     926        25398 :                nunique = nunique + 1
     927        25398 :                iend = j - 1
     928              :                CALL set_potparm_index(potparm, my_index(istart:iend), pot_target, &
     929        25398 :                                       ntype, tmpij, atomic_kind_set, shift_cutoff, do_zbl)
     930              :                !
     931       496115 :                DO i = istart, iend
     932       470717 :                   locij = my_index(i)
     933       470717 :                   CALL get_indexes(locij, ntype, tmpij)
     934       470717 :                   tmp_index(tmpij(1), tmpij(2)) = nunique
     935       496115 :                   tmp_index(tmpij(2), tmpij(1)) = nunique
     936              :                END DO
     937        25398 :                istart = j
     938        25398 :                locij = my_index(j)
     939       489916 :                CALL get_indexes(locij, ntype, tmpij0)
     940              :             END DO
     941         4096 :             nunique = nunique + 1
     942         4096 :             iend = ndim
     943              :             CALL set_potparm_index(potparm, my_index(istart:iend), pot_target, &
     944         4096 :                                    ntype, tmpij, atomic_kind_set, shift_cutoff, do_zbl)
     945        19199 :             DO i = istart, iend
     946        15103 :                locij = my_index(i)
     947        15103 :                CALL get_indexes(locij, ntype, tmpij)
     948        15103 :                tmp_index(tmpij(1), tmpij(2)) = nunique
     949        19199 :                tmp_index(tmpij(2), tmpij(1)) = nunique
     950              :             END DO
     951         4096 :             DEALLOCATE (Cwork)
     952         4096 :             DEALLOCATE (pot_par)
     953              :          ELSE
     954         2063 :             nunique = nunique + 1
     955              :             CALL set_potparm_index(potparm, my_index, pot_target, ntype, tmpij, &
     956         2063 :                                    atomic_kind_set, shift_cutoff, do_zbl)
     957              :          END IF
     958       120842 :          DEALLOCATE (my_index)
     959              :       END DO
     960              :       ! Multiple defined potential
     961              :       n = 0
     962        27664 :       DO i = 1, ntype
     963       543092 :          DO j = 1, i
     964       515428 :             n = n + 1
     965       515428 :             IF (SIZE(potparm%pot(i, j)%pot%type) == 1) CYCLE
     966            8 :             nunique = nunique + 1
     967            8 :             tmp_index(i, j) = nunique
     968            8 :             tmp_index(j, i) = nunique
     969              :             !
     970              :             CALL set_potparm_index(potparm, (/n/), multi_type, ntype, tmpij, &
     971       537846 :                                    atomic_kind_set, shift_cutoff, do_zbl)
     972              :          END DO
     973              :       END DO
     974              :       ! Concluding the postprocess..
     975         5254 :       ALLOCATE (spline_env)
     976         5254 :       CALL spline_env_create(spline_env, ntype, nunique)
     977      1036110 :       spline_env%spltab = tmp_index
     978         5254 :       DEALLOCATE (tmp_index)
     979         5254 :       CALL timestop(handle)
     980        15762 :    END SUBROUTINE get_nonbond_storage
     981              : 
     982              : ! **************************************************************************************************
     983              : !> \brief Trivial for non LJ potential.. gives back in the case of LJ
     984              : !>      the potparm with the smallest sigma..
     985              : !> \param potparm ...
     986              : !> \param my_index ...
     987              : !> \param pot_target ...
     988              : !> \param ntype ...
     989              : !> \param tmpij_out ...
     990              : !> \param atomic_kind_set ...
     991              : !> \param shift_cutoff ...
     992              : !> \param do_zbl ...
     993              : !> \author Teodoro Laino [tlaino] 2007.06
     994              : ! **************************************************************************************************
     995        31565 :    SUBROUTINE set_potparm_index(potparm, my_index, pot_target, ntype, tmpij_out, &
     996              :                                 atomic_kind_set, shift_cutoff, do_zbl)
     997              : 
     998              :       TYPE(pair_potential_pp_type), POINTER              :: potparm
     999              :       INTEGER, INTENT(IN)                                :: my_index(:), pot_target, ntype
    1000              :       INTEGER, INTENT(OUT)                               :: tmpij_out(2)
    1001              :       TYPE(atomic_kind_type), DIMENSION(:), POINTER      :: atomic_kind_set
    1002              :       LOGICAL, INTENT(IN)                                :: shift_cutoff, do_zbl
    1003              : 
    1004              :       CHARACTER(len=*), PARAMETER                        :: routineN = 'set_potparm_index'
    1005              : 
    1006              :       INTEGER                                            :: handle, i, min_val, nvalues, tmpij(2), &
    1007              :                                                             value, zi, zj
    1008        31565 :       INTEGER, ALLOCATABLE, DIMENSION(:)                 :: wrk
    1009              :       LOGICAL                                            :: check
    1010              :       REAL(KIND=dp)                                      :: hicut0, l_epsilon, l_sigma6, m_epsilon, &
    1011              :                                                             m_sigma6, min_sigma6, rcovi, rcovj
    1012        31565 :       REAL(KIND=dp), ALLOCATABLE, DIMENSION(:)           :: sigma6
    1013              :       TYPE(atomic_kind_type), POINTER                    :: atomic_kind
    1014              :       TYPE(pair_potential_single_type), POINTER          :: pot, pot_ref
    1015              : 
    1016        31565 :       CALL timeset(routineN, handle)
    1017              : 
    1018        31565 :       NULLIFY (pot, pot_ref)
    1019        31565 :       nvalues = SIZE(my_index)
    1020        31565 :       IF ((pot_target == lj_type) .OR. (pot_target == lj_charmm_type)) THEN
    1021        74913 :          ALLOCATE (sigma6(nvalues))
    1022        74913 :          ALLOCATE (wrk(nvalues))
    1023       505927 :          min_sigma6 = HUGE(0.0_dp)
    1024       505927 :          m_epsilon = -HUGE(0.0_dp)
    1025       505927 :          DO i = 1, nvalues
    1026       480956 :             value = my_index(i)
    1027       480956 :             CALL get_indexes(value, ntype, tmpij)
    1028       480956 :             pot => potparm%pot(tmpij(1), tmpij(2))%pot
    1029              :             ! Preliminary check..
    1030       480956 :             check = SIZE(pot%type) == 1
    1031       480956 :             CPASSERT(check)
    1032              : 
    1033       480956 :             sigma6(i) = pot%set(1)%lj%sigma6
    1034       480956 :             l_epsilon = pot%set(1)%lj%epsilon
    1035       480956 :             IF (sigma6(i) /= 0.0_dp) min_sigma6 = MIN(min_sigma6, sigma6(i))
    1036       480956 :             IF (sigma6(i) == 0.0_dp) sigma6(i) = -HUGE(0.0_dp)
    1037       505927 :             IF (l_epsilon /= 0.0_dp) m_epsilon = MAX(m_epsilon, l_epsilon)
    1038              :          END DO
    1039        24971 :          CALL sort(sigma6, nvalues, wrk)
    1040        24971 :          min_val = my_index(wrk(nvalues))
    1041        24971 :          m_sigma6 = sigma6(nvalues)
    1042              :          ! In case there are only zeros.. let's consider them properly..
    1043        24971 :          IF (m_sigma6 == -HUGE(0.0_dp)) m_sigma6 = 1.0_dp
    1044        24971 :          IF (m_epsilon == -HUGE(0.0_dp)) m_epsilon = 0.0_dp
    1045        24971 :          IF (min_sigma6 == HUGE(0.0_dp)) min_sigma6 = 0.0_dp
    1046        24971 :          DEALLOCATE (sigma6)
    1047        24971 :          DEALLOCATE (wrk)
    1048              :       ELSE
    1049        41066 :          min_val = MINVAL(my_index(:))
    1050              :       END IF
    1051        31565 :       CALL get_indexes(min_val, ntype, tmpij)
    1052        31565 :       tmpij_out = tmpij
    1053        31565 :       pot => potparm%pot(tmpij(1), tmpij(2))%pot
    1054        31565 :       pot%undef = .TRUE.
    1055        31565 :       IF (shift_cutoff) THEN
    1056        28045 :          hicut0 = SQRT(pot%rcutsq)
    1057        28045 :          IF (ABS(hicut0) <= MIN_HICUT_VALUE) hicut0 = DEFAULT_HICUT_VALUE
    1058              :       END IF
    1059        31565 :       CALL init_genpot(potparm, ntype)
    1060              : 
    1061       546993 :       DO i = 1, nvalues
    1062       515428 :          value = my_index(i)
    1063       515428 :          CALL get_indexes(value, ntype, tmpij)
    1064       515428 :          pot => potparm%pot(tmpij(1), tmpij(2))%pot
    1065       515428 :          CALL spline_factor_create(pot%spl_f)
    1066       515428 :          pot%spl_f%rcutsq_f = 1.0_dp
    1067      1030856 :          pot%spl_f%rscale = 1.0_dp
    1068      1062421 :          pot%spl_f%fscale = 1.0_dp
    1069              :       END DO
    1070              : 
    1071        94691 :       IF (ANY(potential_single_allocation == pot_target)) THEN
    1072         9388 :          DO i = 1, nvalues
    1073         9384 :             value = my_index(i)
    1074         9384 :             CALL get_indexes(value, ntype, tmpij)
    1075         9384 :             pot => potparm%pot(tmpij(1), tmpij(2))%pot
    1076              : 
    1077         9384 :             check = SIZE(pot%type) == 1
    1078         9384 :             CPASSERT(check)
    1079              :             ! Undef potential.. this will be used to compute the splines..
    1080         9388 :             IF ((pot_target == lj_type) .OR. (pot_target == lj_charmm_type)) THEN
    1081         9384 :                l_sigma6 = pot%set(1)%lj%sigma6
    1082         9384 :                l_epsilon = pot%set(1)%lj%epsilon
    1083              :                ! Undef potential.. this will be used to compute the splines..
    1084         9384 :                IF (pot%undef) THEN
    1085            4 :                   pot%set(1)%lj%sigma6 = m_sigma6
    1086            4 :                   pot%set(1)%lj%sigma12 = m_sigma6**2
    1087            4 :                   pot%set(1)%lj%epsilon = m_epsilon
    1088              :                END IF
    1089         9384 :                pot%spl_f%rscale(1) = 1.0_dp
    1090         9384 :                pot%spl_f%fscale(1) = 0.0_dp
    1091         9384 :                IF (l_sigma6*l_epsilon /= 0.0_dp) THEN
    1092         9384 :                   pot%spl_f%rcutsq_f = (min_sigma6/m_sigma6)**(1.0_dp/3.0_dp)
    1093         9384 :                   pot%spl_f%rscale(1) = (l_sigma6/m_sigma6)**(1.0_dp/3.0_dp)
    1094         9384 :                   pot%spl_f%fscale(1) = l_epsilon/m_epsilon
    1095              :                END IF
    1096              :             END IF
    1097              :          END DO
    1098              :       END IF
    1099              : 
    1100       546993 :       DO i = 1, nvalues
    1101       515428 :          value = my_index(i)
    1102       515428 :          CALL get_indexes(value, ntype, tmpij)
    1103       515428 :          pot => potparm%pot(tmpij(1), tmpij(2))%pot
    1104              : 
    1105       515428 :          IF (do_zbl) THEN
    1106           48 :             atomic_kind => atomic_kind_set(tmpij(1))
    1107           48 :             CALL get_atomic_kind(atomic_kind, rcov=rcovi, z=zi)
    1108           48 :             atomic_kind => atomic_kind_set(tmpij(2))
    1109           48 :             CALL get_atomic_kind(atomic_kind, rcov=rcovj, z=zj)
    1110              :             CALL zbl_matching_polinomial(pot, rcovi, rcovj, REAL(zi, KIND=dp), &
    1111           48 :                                          REAL(zj, KIND=dp))
    1112              :          END IF
    1113              :          ! Derivative factors
    1114      1030856 :          pot%spl_f%dscale = pot%spl_f%fscale/pot%spl_f%rscale
    1115              :          ! Cutoff for the potentials on splines
    1116       546993 :          IF (shift_cutoff) THEN
    1117              :             ! Cutoff NonBonded
    1118       115008 :             pot%spl_f%cutoff = ener_pot(pot, hicut0, 0.0_dp)
    1119              :          END IF
    1120              :       END DO
    1121              : 
    1122              :       ! Handle the cutoff
    1123        31565 :       IF (shift_cutoff) THEN
    1124        28045 :          pot_ref => potparm%pot(tmpij_out(1), tmpij_out(2))%pot
    1125       143053 :          DO i = 1, nvalues
    1126       115008 :             value = my_index(i)
    1127       115008 :             CALL get_indexes(value, ntype, tmpij)
    1128       115008 :             pot => potparm%pot(tmpij(1), tmpij(2))%pot
    1129       115008 :             IF (value == min_val) CYCLE
    1130              :             ! Cutoff NonBonded
    1131       143053 :             pot%spl_f%cutoff = pot_ref%spl_f%cutoff*pot%spl_f%fscale(1) - pot%spl_f%cutoff
    1132              :          END DO
    1133              :       END IF
    1134        31565 :       CALL finalizef()
    1135              : 
    1136        31565 :       CALL timestop(handle)
    1137              : 
    1138        31565 :    END SUBROUTINE set_potparm_index
    1139              : 
    1140              : ! **************************************************************************************************
    1141              : !> \brief Gives back the indices of the matrix w.r.t. the collective array index
    1142              : !> \param Inind ...
    1143              : !> \param ndim ...
    1144              : !> \param ij ...
    1145              : !> \author Teodoro Laino [tlaino] 2006.05
    1146              : ! **************************************************************************************************
    1147      2668903 :    SUBROUTINE get_indexes(Inind, ndim, ij)
    1148              :       INTEGER, INTENT(IN)                                :: Inind, ndim
    1149              :       INTEGER, DIMENSION(2), INTENT(OUT)                 :: ij
    1150              : 
    1151              :       INTEGER                                            :: i, tmp
    1152              : 
    1153      2668903 :       tmp = 0
    1154      8006709 :       ij = HUGE(0)
    1155    355422195 :       DO i = 1, ndim
    1156    355422195 :          tmp = tmp + i
    1157    355422195 :          IF (tmp >= Inind) THEN
    1158      2668903 :             ij(1) = i
    1159      2668903 :             ij(2) = Inind - tmp + i
    1160      2668903 :             EXIT
    1161              :          END IF
    1162              :       END DO
    1163      2668903 :    END SUBROUTINE get_indexes
    1164              : 
    1165              : END MODULE pair_potential
    1166              : 
        

Generated by: LCOV version 2.0-1