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