(git:e5b1968)
Loading...
Searching...
No Matches
pair_potential.F
Go to the documentation of this file.
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! **************************************************************************************************
16MODULE 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
69CONTAINS
70
71! **************************************************************************************************
72!> \brief Initialize genpot
73!> \param potparm ...
74!> \param ntype ...
75!> \par History
76!> Teo 2007.06 - Zurich University
77! **************************************************************************************************
78 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 CALL timeset(routinen, handle)
88
89 NULLIFY (pot)
90 ngp = 0
91 ! Prescreen for general potential type
92 DO i = 1, ntype ! i: first atom type
93 DO j = 1, i ! j: second atom type
94 pot => potparm%pot(i, j)%pot
95 ngp = ngp + count(pot%type == gp_type)
96 END DO
97 END DO
98 CALL initf(ngp)
99 ngp = 0
100 DO i = 1, ntype ! i: first atom type
101 DO j = 1, i ! j: second atom type
102 pot => potparm%pot(i, j)%pot
103 DO k = 1, SIZE(pot%type)
104 IF (pot%type(k) == gp_type) THEN
105 ngp = ngp + 1
106 pot%set(k)%gp%myid = ngp
107 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 CALL timestop(handle)
113
114 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 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 n = 0
157 ncount = 0
158
159 ntype = SIZE(atomic_kind_set)
160 CALL timeset(routinen, handle)
161 IF (iw3 > 0) THEN
162 WRITE (iw3, "(/,T2,A,I0,A,I0,A)") &
163 "SPLINE_INFO| Generating ", (ntype*(ntype + 1))/2, " splines for "// &
164 trim(adjustl(nonbonded_type))//" interactions "
165 WRITE (iw3, "(T2,A,I0,A)") &
166 " Due to ", ntype, " different atomic kinds"
167 END IF
168 CALL init_genpot(potparm, ntype)
169 ! Real computation of splines
170 ip = 0
171 DO i = 1, ntype
172 DO j = 1, i
173 pot => potparm%pot(i, j)%pot
174 IF (iw3 > 0 .AND. iw <= 0) THEN
175 IF (mod(i*(i - 1)/2 + j, max(1, (ntype*(ntype + 1))/(2*10))) == 0) THEN
176 WRITE (unit=iw3, advance="NO", fmt='(2X,A3,I0)') '...', i*(i - 1)/2 + j
177 ip = ip + 1
178 IF (ip >= 11) THEN
179 WRITE (iw3, *)
180 ip = 0
181 END IF
182 END IF
183 END IF
184 ! Setup of Exclusion Types
185 pot%no_pp = .true.
186 pot%no_mb = .true.
187 DO k = 1, SIZE(pot%type)
188 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 pot%no_pp = .false.
192 CASE (tersoff_type)
193 pot%no_mb = .false.
194 CASE (siepmann_type)
195 pot%no_mb = .false.
196 CASE (gal_type)
197 pot%no_mb = .false.
198 CASE (gal21_type)
199 pot%no_mb = .false.
200 CASE (nn_type)
201 ! Do nothing..
202 CASE DEFAULT
203 ! Never reach this point
204 cpabort("")
205 END SELECT
206 ! Special case for EAM
207 SELECT CASE (pot%type(k))
208 CASE (ea_type, quip_type, nequip_type, allegro_type, deepmd_type, ace_type)
209 pot%no_mb = .false.
210 END SELECT
211 END DO
212
213 ! Starting SetUp of splines
214 IF (.NOT. pot%undef) cycle
215 ncount = ncount + 1
216 n = spline_env%spltab(i, j)
217 locut = rlow_nb
218 hicut0 = sqrt(pot%rcutsq)
219 IF (abs(hicut0) <= min_hicut_value) hicut0 = default_hicut_value
220 hicut = hicut0/sqrt(pot%spl_f%rcutsq_f)
221
222 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 energy_cutoff, emax_spline)
227 locut = max(locut*sqrt(pot%spl_f%rcutsq_f), rlow_nb)
228
229 ! Real Generation of the Spline
230 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 nonbonded_type)
235
236 pot%undef = .false.
237 ! Unique Spline working only for a pure LJ potential..
238 IF (SIZE(pot%type) == 1) THEN
239 IF (any(potential_single_allocation == pot%type(1))) THEN
240 ! Restoring the proper values for the generating spline pot
241 IF ((pot%type(1) == lj_type) .OR. (pot%type(1) == lj_charmm_type)) THEN
242 pot%set(1)%lj%sigma6 = pot%set(1)%lj%sigma6*pot%spl_f%rscale(1)**3
243 pot%set(1)%lj%sigma12 = pot%set(1)%lj%sigma6**2
244 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 IF (shift_cutoff) THEN
250 pot%spl_f%cutoff = pot%spl_f%cutoff*pot%spl_f%fscale(1) - &
251 ener_pot(pot, hicut0, 0.0_dp)
252 END IF
253 END DO
254 END DO
255 CALL finalizef()
256
257 IF (iw > 0) THEN
258 WRITE (iw, '(/,T2,A,I0)') &
259 "SPLINE_INFO| Number of pair potential splines allocated: ", maxval(spline_env%spltab)
260 END IF
261 IF (iw3 > 0) THEN
262 WRITE (iw3, '(/,T2,A,I0,/,T2,A)') &
263 "SPLINE_INFO| Number of unique splines computed: ", maxval(spline_env%spltab), &
264 "SPLINE_INFO| Done"
265 END IF
266
267 CALL timestop(handle)
268
269 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 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 dx2 = (hicut - locut)
300 x = hicut
301 locut_found = locut
302 found_locut = .false.
303 DO ilevel = 1, 2
304 dx2 = dx2/100.0_dp
305 DO jx = 1, 100
306 e = ener_pot(pot, x, energy_cutoff)
307 IF (do_zbl) THEN
308 e = e + ener_zbl(pot, x)
309 END IF
310 IF (abs(e) > emax_spline) THEN
311 locut_found = x
312 found_locut = .true.
313 EXIT
314 END IF
315 x = x - dx2
316 END DO
317 x = x + dx2
318 END DO
319 locut = locut_found
320
321 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 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 NULLIFY (logger, spl_f)
374 logger => cp_get_default_logger()
375
376 CALL spline_factor_create(spl_f)
377 mfac = 5
378 IF (npoints > 0) THEN
379 fixed_spline_points = .true.
380 ELSE
381 fixed_spline_points = .false.
382 npoints = 20
383 IF (.NOT. found_locut) npoints = 2
384 END IF
385 spline_data => spl_p(1)%spline_data
386 DO WHILE (.true.)
387 CALL init_splinexy(spline_data, npoints + 1)
388 dx2 = (1.0_dp/locut**2 - 1.0_dp/hicut**2)/real(npoints, kind=dp)
389 x2 = 1.0_dp/hicut**2
390 spline_data%x1 = x2
391 DO jx = 1, npoints + 1
392 ! jx: loop over 1/distance**2
393 x = sqrt(1.0_dp/x2)
394 e = ener_pot(pot, x, energy_cutoff)
395 IF (do_zbl) THEN
396 e = e + ener_zbl(pot, x)
397 END IF
398 spline_data%y(jx) = e
399 x2 = x2 + dx2
400 END DO
401 CALL init_spline(spline_data, dx=dx2)
402 ! This is the check for required accuracy on spline setup
403 dx2 = (hicut - locut)/real(mfac*npoints + 1, kind=dp)
404 x2 = locut + dx2
405 diffmax = -1.0_dp
406 xsav = hicut
407 ! if a fixed number of points is requested, no check on its error
408 IF (fixed_spline_points) EXIT
409 DO jx = 1, mfac*npoints
410 x = x2
411 e = ener_pot(pot, x, energy_cutoff)
412 IF (do_zbl) THEN
413 e = e + ener_zbl(pot, x)
414 END IF
415 IF (abs(e) < max_energy) THEN
416 xdum1 = abs(e - potential_s(spl_p, x*x, xdum, spl_f, logger))
417 diffmax = max(diffmax, xdum1)
418 xsav = min(x, xsav)
419 END IF
420 x2 = x2 + dx2
421 IF (x2 > hicut) EXIT
422 END DO
423 IF (npoints > max_points) THEN
424 WRITE (message, '(A,I8,A,G12.6,A)') "SPLINE_INFO| Number of points: ", npoints, &
425 " obtained accuracy ", diffmax, ". MM SPLINE: no convergence on required"// &
426 " accuracy (adjust EPS_SPLINE and rerun)"
427 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 IF (diffmax > eps_spline .OR. diffmax < 0.0_dp) THEN
431 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 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 " SPLINE_INFO| Spline number: ", ncount, &
441 " SPLINE_INFO| Unique spline number: ", n, &
442 " SPLINE_INFO| Atomic kind numbers: ", i, j, &
443 " SPLINE_INFO| Atomic kind names: "//trim(adjustl(atomic_kind_set(i)%name))//" "// &
444 trim(adjustl(atomic_kind_set(j)%name)), &
445 " SPLINE_INFO| Number of spline points: ", npoints, &
446 " SPLINE_INFO| Requested accuracy [Hartree]: ", eps_spline, &
447 " SPLINE_INFO| Achieved accuracy [Hartree]: ", diffmax, &
448 " SPLINE_INFO| Spline range [bohr]: ", locut, hicut, &
449 " SPLINE_INFO| Spline range used to achieve accuracy [bohr]:", xsav, hicut
450 dx2 = (hicut - locut)/real(npoints + 1, kind=dp)
451 x = locut + dx2
452 WRITE (unit=iw, fmt='(A,ES17.9)') &
453 " SPLINE_INFO| Spline value at RMIN [Hartree]: ", potential_s(spl_p, x*x, xdum, spl_f, logger), &
454 " SPLINE_INFO| Spline value at RMAX [Hartree]: ", potential_s(spl_p, hicut*hicut, xdum, spl_f, logger), &
455 " SPLINE_INFO| Non-bonded energy cutoff [Hartree]: ", energy_cutoff
456 END IF
457 ! Print spline data on file if requested
458 IF (iw2 > 0) THEN
459 ! Set increment to 200 points per Angstrom
460 nppa = 200
461 dx = bohr/real(nppa, kind=dp)
462 nx = nint(hicut/dx)
463 file_name = ""
464 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 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 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 "# Spline number: ", ncount, &
477 "# Unique spline number: ", n, &
478 "# Atomic kind numbers: ", i, j, &
479 "# Atomic kind names: "//trim(adjustl(atomic_kind_set(i)%name))//" "// &
480 trim(adjustl(atomic_kind_set(j)%name)), &
481 "# Number of spline points: ", npoints, &
482 "# Requested accuracy [eV]: ", eps_spline*evolt, &
483 "# Achieved accuracy [eV]: ", diffmax*evolt, &
484 "# Spline range [Angstrom]: ", locut/bohr, hicut/bohr, &
485 "# Spline range used to achieve accuracy [Angstrom]:", xsav/bohr, hicut/bohr, &
486 "# Non-bonded energy cutoff [eV]: ", energy_cutoff*evolt, &
487 "# Test spline using ", nppa, " points per Angstrom:", &
488 "# Abscissa [Angstrom] Energy [eV] Splined energy [eV] Derivative [eV/Angstrom]"// &
489 " |Energy error| [eV]"
490 x = 0.0_dp
491 DO jx = 0, nx
492 IF (x > hicut) x = hicut
493 IF (x > locut) THEN
494 e = ener_pot(pot, x, energy_cutoff)
495 IF (do_zbl) e = e + ener_zbl(pot, x)
496 e_spline = potential_s(spl_p, x*x, xdum, spl_f, logger)
497 WRITE (unit=unit_number, fmt="(5ES25.12)") &
498 x/bohr, e*evolt, e_spline*evolt, -bohr*x*xdum*evolt, abs((e - e_spline)*evolt)
499 END IF
500 x = x + dx
501 END DO
502 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 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 unit_number=unit_number)
513 ! Recommended increment for dp is 0.0005 nm = 0.005 Angstrom
514 ! which are 200 points/Angstrom
515 rcut = 0.1_dp*hicut/bohr
516 x = 0.0_dp
517 DO jx = 0, nx
518 IF (x > hicut) x = hicut
519 r = 0.1_dp*x/bohr ! Convert bohr to nm
520 IF (x <= locut) THEN
521 WRITE (unit=unit_number, fmt="(7ES25.12)") &
522 r, (0.0_dp, ix=1, 6)
523 ELSE
524 e_spline = potential_s(spl_p, x*x, xdum, spl_f, logger)
525 f = 1.0_dp/r
526 df = -1.0_dp/r**2
527 g = -1.0_dp/r**6 + 1.0_dp/rcut**6
528 dg = 6.0_dp/r**7
529 h = e_spline*kjmol
530 dh = -10.0_dp*bohr*x*xdum*kjmol
531 WRITE (unit=unit_number, fmt="(7ES25.12)") &
532 r, f, -df, & ! r, f(r), -f'(r) => probably not used
533 g, -dg, & ! g(r), -g'(r) => not used, if C = 0
534 h, -dh ! h(r), -h'(r) => used, if A = 1
535 END IF
536 x = x + dx
537 END DO
538 CALL close_file(unit_number=unit_number)
539 END IF
540
541 CALL spline_factor_release(spl_f)
542
543 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 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 INTEGER, ALLOCATABLE, DIMENSION(:) :: iwork1, iwork2, my_index
568 INTEGER, ALLOCATABLE, DIMENSION(:, :) :: tmp_index
569 LOGICAL :: at_least_one, check
570 REAL(kind=dp), ALLOCATABLE, DIMENSION(:) :: cwork, rwork, wtmp
571 REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :) :: pot_par
572
573 CALL timeset(routinen, handle)
574
575 ntype = SIZE(atomic_kind_set)
576 DO i = 1, ntype
577 DO j = 1, i
578 potparm%pot(i, j)%pot%undef = .false.
579 END DO
580 END DO
581 ALLOCATE (tmp_index(ntype, ntype))
582 !
583 nunique = 0
584 tmp_index = huge(0)
585 DO pot_target = minval(list_pot), maxval(list_pot)
586 ndim = 0
587 DO i = 1, ntype
588 DO j = 1, i
589 IF (SIZE(potparm%pot(i, j)%pot%type) /= 1) cycle
590 IF (potparm%pot(i, j)%pot%type(1) == pot_target) THEN
591 tmp_index(i, j) = 1
592 tmp_index(j, i) = 1
593 ndim = ndim + 1
594 END IF
595 END DO
596 END DO
597 IF (ndim == 0) cycle ! No potential of this kind found
598 nvar = 0
599 SELECT CASE (pot_target)
600 CASE (lj_type, lj_charmm_type)
601 nvar = 3 + nvar
602 CASE (wl_type)
603 nvar = 3 + nvar
604 CASE (gw_type)
605 nvar = 5 + nvar
606 CASE (ea_type)
607 nvar = 4 + nvar
608 CASE (quip_type)
609 nvar = 1 + nvar
610 CASE (nequip_type)
611 nvar = 1 + nvar
612 CASE (allegro_type)
613 nvar = 1 + nvar
614 CASE (ace_type)
615 nvar = 2 + nvar
616 CASE (deepmd_type)
617 nvar = 2 + nvar
618 CASE (ft_type)
619 nvar = 4 + nvar
620 CASE (ftd_type)
621 nvar = 6 + nvar
622 CASE (ip_type)
623 nvar = 3 + nvar
624 CASE (b4_type)
625 nvar = 6 + nvar
626 CASE (bm_type)
627 nvar = 9 + nvar
628 CASE (gp_type)
629 nvar = 2 + nvar
630 CASE (tersoff_type)
631 nvar = 13 + nvar
632 CASE (siepmann_type)
633 nvar = 5 + nvar
634 CASE (gal_type)
635 nvar = 12 + nvar
636 CASE (gal21_type)
637 nvar = 30 + nvar
638 CASE (nn_type)
639 nvar = nvar
640 CASE (tab_type)
641 nvar = 4 + nvar
642 CASE DEFAULT
643 cpabort("")
644 END SELECT
645 ! Setup a table of the indexes..
646 ALLOCATE (my_index(ndim))
647 n = 0
648 nk = 0
649 DO i = 1, ntype
650 DO j = 1, i
651 n = n + 1
652 IF (SIZE(potparm%pot(i, j)%pot%type) /= 1) cycle
653 IF (potparm%pot(i, j)%pot%type(1) == pot_target) THEN
654 nk = nk + 1
655 my_index(nk) = n
656 END IF
657 END DO
658 END DO
659 IF (nvar /= 0) THEN
660 ALLOCATE (pot_par(ndim, nvar))
661 n = 0
662 nk = 0
663 DO i = 1, ntype
664 DO j = 1, i
665 n = n + 1
666 IF (SIZE(potparm%pot(i, j)%pot%type) /= 1) cycle
667 IF (potparm%pot(i, j)%pot%type(1) == pot_target) THEN
668 nk = nk + 1
669 my_index(nk) = n
670 SELECT CASE (pot_target)
671 CASE (lj_type, lj_charmm_type)
672 pot_par(nk, 1) = potparm%pot(i, j)%pot%set(1)%lj%epsilon
673 pot_par(nk, 2) = potparm%pot(i, j)%pot%set(1)%lj%sigma6
674 pot_par(nk, 3) = potparm%pot(i, j)%pot%set(1)%lj%sigma12
675 CASE (gp_type)
676 pot_par(nk, 1) = str2id(potparm%pot(i, j)%pot%set(1)%gp%potential)
677 pot_par(nk, 2) = str2id(potparm%pot(i, j)%pot%set(1)%gp%variables)
678 CASE (wl_type)
679 pot_par(nk, 1) = potparm%pot(i, j)%pot%set(1)%willis%a
680 pot_par(nk, 2) = potparm%pot(i, j)%pot%set(1)%willis%b
681 pot_par(nk, 3) = potparm%pot(i, j)%pot%set(1)%willis%c
682 CASE (gw_type)
683 pot_par(nk, 1) = potparm%pot(i, j)%pot%set(1)%goodwin%vr0
684 pot_par(nk, 2) = potparm%pot(i, j)%pot%set(1)%goodwin%m
685 pot_par(nk, 3) = potparm%pot(i, j)%pot%set(1)%goodwin%mc
686 pot_par(nk, 4) = potparm%pot(i, j)%pot%set(1)%goodwin%d
687 pot_par(nk, 5) = potparm%pot(i, j)%pot%set(1)%goodwin%dc
688 CASE (ea_type)
689 pot_par(nk, 1) = potparm%pot(i, j)%pot%set(1)%eam%drar
690 pot_par(nk, 2) = potparm%pot(i, j)%pot%set(1)%eam%drhoar
691 pot_par(nk, 3) = potparm%pot(i, j)%pot%set(1)%eam%acutal
692 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 trim(potparm%pot(i, j)%pot%set(1)%quip%calc_args))
698 CASE (nequip_type)
699 pot_par(nk, 1) = str2id( &
700 trim(potparm%pot(i, j)%pot%set(1)%nequip%nequip_file_name))
701 CASE (allegro_type)
702 pot_par(nk, 1) = str2id( &
703 trim(potparm%pot(i, j)%pot%set(1)%allegro%allegro_file_name))
704 CASE (ace_type)
705 pot_par(nk, 1) = str2id( &
706 trim(potparm%pot(i, j)%pot%set(1)%ace%ace_file_name))
707 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 trim(potparm%pot(i, j)%pot%set(1)%deepmd%deepmd_file_name))
711 pot_par(nk, 2) = potparm%pot(i, j)%pot%set(1)%deepmd%atom_deepmd_type
712 CASE (ft_type)
713 pot_par(nk, 1) = potparm%pot(i, j)%pot%set(1)%ft%A
714 pot_par(nk, 2) = potparm%pot(i, j)%pot%set(1)%ft%B
715 pot_par(nk, 3) = potparm%pot(i, j)%pot%set(1)%ft%C
716 pot_par(nk, 4) = potparm%pot(i, j)%pot%set(1)%ft%D
717 CASE (ftd_type)
718 pot_par(nk, 1) = potparm%pot(i, j)%pot%set(1)%ftd%A
719 pot_par(nk, 2) = potparm%pot(i, j)%pot%set(1)%ftd%B
720 pot_par(nk, 3) = potparm%pot(i, j)%pot%set(1)%ftd%C
721 pot_par(nk, 4) = potparm%pot(i, j)%pot%set(1)%ftd%D
722 pot_par(nk, 5) = potparm%pot(i, j)%pot%set(1)%ftd%BD(1)
723 pot_par(nk, 6) = potparm%pot(i, j)%pot%set(1)%ftd%BD(2)
724 CASE (ip_type)
725 pot_par(nk, 1) = potparm%pot(i, j)%pot%set(1)%ipbv%rcore
726 pot_par(nk, 2) = potparm%pot(i, j)%pot%set(1)%ipbv%m
727 pot_par(nk, 3) = potparm%pot(i, j)%pot%set(1)%ipbv%b
728 CASE (b4_type)
729 pot_par(nk, 1) = potparm%pot(i, j)%pot%set(1)%buck4r%a
730 pot_par(nk, 2) = potparm%pot(i, j)%pot%set(1)%buck4r%b
731 pot_par(nk, 3) = potparm%pot(i, j)%pot%set(1)%buck4r%c
732 pot_par(nk, 4) = potparm%pot(i, j)%pot%set(1)%buck4r%r1
733 pot_par(nk, 5) = potparm%pot(i, j)%pot%set(1)%buck4r%r2
734 pot_par(nk, 6) = potparm%pot(i, j)%pot%set(1)%buck4r%r3
735 CASE (bm_type)
736 pot_par(nk, 1) = potparm%pot(i, j)%pot%set(1)%buckmo%f0
737 pot_par(nk, 2) = potparm%pot(i, j)%pot%set(1)%buckmo%a1
738 pot_par(nk, 3) = potparm%pot(i, j)%pot%set(1)%buckmo%a2
739 pot_par(nk, 4) = potparm%pot(i, j)%pot%set(1)%buckmo%b1
740 pot_par(nk, 5) = potparm%pot(i, j)%pot%set(1)%buckmo%b2
741 pot_par(nk, 6) = potparm%pot(i, j)%pot%set(1)%buckmo%c
742 pot_par(nk, 7) = potparm%pot(i, j)%pot%set(1)%buckmo%d
743 pot_par(nk, 8) = potparm%pot(i, j)%pot%set(1)%buckmo%r0
744 pot_par(nk, 9) = potparm%pot(i, j)%pot%set(1)%buckmo%beta
745 CASE (tersoff_type)
746 pot_par(nk, 1) = potparm%pot(i, j)%pot%set(1)%tersoff%A
747 pot_par(nk, 2) = potparm%pot(i, j)%pot%set(1)%tersoff%B
748 pot_par(nk, 3) = potparm%pot(i, j)%pot%set(1)%tersoff%lambda1
749 pot_par(nk, 4) = potparm%pot(i, j)%pot%set(1)%tersoff%lambda2
750 pot_par(nk, 5) = potparm%pot(i, j)%pot%set(1)%tersoff%alpha
751 pot_par(nk, 6) = potparm%pot(i, j)%pot%set(1)%tersoff%beta
752 pot_par(nk, 7) = potparm%pot(i, j)%pot%set(1)%tersoff%n
753 pot_par(nk, 8) = potparm%pot(i, j)%pot%set(1)%tersoff%c
754 pot_par(nk, 9) = potparm%pot(i, j)%pot%set(1)%tersoff%d
755 pot_par(nk, 10) = potparm%pot(i, j)%pot%set(1)%tersoff%h
756 pot_par(nk, 11) = potparm%pot(i, j)%pot%set(1)%tersoff%lambda3
757 pot_par(nk, 12) = potparm%pot(i, j)%pot%set(1)%tersoff%bigR
758 pot_par(nk, 13) = potparm%pot(i, j)%pot%set(1)%tersoff%bigD
759 CASE (siepmann_type)
760 pot_par(nk, 1) = potparm%pot(i, j)%pot%set(1)%siepmann%B
761 pot_par(nk, 2) = potparm%pot(i, j)%pot%set(1)%siepmann%D
762 pot_par(nk, 3) = potparm%pot(i, j)%pot%set(1)%siepmann%E
763 pot_par(nk, 4) = potparm%pot(i, j)%pot%set(1)%siepmann%F
764 pot_par(nk, 5) = potparm%pot(i, j)%pot%set(1)%siepmann%beta
765 CASE (gal_type)
766 pot_par(nk, 1) = potparm%pot(i, j)%pot%set(1)%gal%epsilon
767 pot_par(nk, 2) = potparm%pot(i, j)%pot%set(1)%gal%bxy
768 pot_par(nk, 3) = potparm%pot(i, j)%pot%set(1)%gal%bz
769 pot_par(nk, 4) = potparm%pot(i, j)%pot%set(1)%gal%r1
770 pot_par(nk, 5) = potparm%pot(i, j)%pot%set(1)%gal%r2
771 pot_par(nk, 6) = potparm%pot(i, j)%pot%set(1)%gal%a1
772 pot_par(nk, 7) = potparm%pot(i, j)%pot%set(1)%gal%a2
773 pot_par(nk, 8) = potparm%pot(i, j)%pot%set(1)%gal%a3
774 pot_par(nk, 9) = potparm%pot(i, j)%pot%set(1)%gal%a4
775 pot_par(nk, 10) = potparm%pot(i, j)%pot%set(1)%gal%a
776 pot_par(nk, 11) = potparm%pot(i, j)%pot%set(1)%gal%b
777 pot_par(nk, 12) = potparm%pot(i, j)%pot%set(1)%gal%c
778 CASE (gal21_type)
779 pot_par(nk, 1) = potparm%pot(i, j)%pot%set(1)%gal21%epsilon1
780 pot_par(nk, 2) = potparm%pot(i, j)%pot%set(1)%gal21%epsilon2
781 pot_par(nk, 3) = potparm%pot(i, j)%pot%set(1)%gal21%epsilon3
782 pot_par(nk, 4) = potparm%pot(i, j)%pot%set(1)%gal21%bxy1
783 pot_par(nk, 5) = potparm%pot(i, j)%pot%set(1)%gal21%bxy2
784 pot_par(nk, 6) = potparm%pot(i, j)%pot%set(1)%gal21%bz1
785 pot_par(nk, 7) = potparm%pot(i, j)%pot%set(1)%gal21%bz2
786 pot_par(nk, 8) = potparm%pot(i, j)%pot%set(1)%gal21%r1
787 pot_par(nk, 9) = potparm%pot(i, j)%pot%set(1)%gal21%r2
788 pot_par(nk, 10) = potparm%pot(i, j)%pot%set(1)%gal21%a11
789 pot_par(nk, 11) = potparm%pot(i, j)%pot%set(1)%gal21%a12
790 pot_par(nk, 12) = potparm%pot(i, j)%pot%set(1)%gal21%a13
791 pot_par(nk, 13) = potparm%pot(i, j)%pot%set(1)%gal21%a21
792 pot_par(nk, 14) = potparm%pot(i, j)%pot%set(1)%gal21%a22
793 pot_par(nk, 15) = potparm%pot(i, j)%pot%set(1)%gal21%a23
794 pot_par(nk, 16) = potparm%pot(i, j)%pot%set(1)%gal21%a31
795 pot_par(nk, 17) = potparm%pot(i, j)%pot%set(1)%gal21%a32
796 pot_par(nk, 18) = potparm%pot(i, j)%pot%set(1)%gal21%a33
797 pot_par(nk, 19) = potparm%pot(i, j)%pot%set(1)%gal21%a41
798 pot_par(nk, 20) = potparm%pot(i, j)%pot%set(1)%gal21%a42
799 pot_par(nk, 21) = potparm%pot(i, j)%pot%set(1)%gal21%a43
800 pot_par(nk, 22) = potparm%pot(i, j)%pot%set(1)%gal21%AO1
801 pot_par(nk, 23) = potparm%pot(i, j)%pot%set(1)%gal21%AO2
802 pot_par(nk, 24) = potparm%pot(i, j)%pot%set(1)%gal21%BO1
803 pot_par(nk, 25) = potparm%pot(i, j)%pot%set(1)%gal21%BO2
804 pot_par(nk, 26) = potparm%pot(i, j)%pot%set(1)%gal21%c
805 pot_par(nk, 27) = potparm%pot(i, j)%pot%set(1)%gal21%AH1
806 pot_par(nk, 28) = potparm%pot(i, j)%pot%set(1)%gal21%AH2
807 pot_par(nk, 29) = potparm%pot(i, j)%pot%set(1)%gal21%BH1
808 pot_par(nk, 30) = potparm%pot(i, j)%pot%set(1)%gal21%BH2
809 CASE (tab_type)
810 pot_par(nk, 1) = potparm%pot(i, j)%pot%set(1)%tab%dr
811 pot_par(nk, 2) = potparm%pot(i, j)%pot%set(1)%tab%rcut
812 pot_par(nk, 3) = potparm%pot(i, j)%pot%set(1)%tab%npoints
813 pot_par(nk, 4) = potparm%pot(i, j)%pot%set(1)%tab%index
814 CASE (nn_type)
815 ! no checks
816 CASE DEFAULT
817 cpabort("")
818 END SELECT
819 IF (any(potential_single_allocation == pot_target)) THEN
820 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 ALLOCATE (rwork(ndim))
827 ALLOCATE (iwork1(ndim))
828 ALLOCATE (iwork2(ndim))
829 ALLOCATE (wtmp(nvar))
830 CALL sort(pot_par(:, 1), ndim, iwork1)
831 ! Sort all the other components of the potential
832 DO k = 2, nvar
833 rwork(:) = pot_par(:, k)
834 DO i = 1, ndim
835 pot_par(i, k) = rwork(iwork1(i))
836 END DO
837 END DO
838 iwork2(:) = my_index
839 DO i = 1, ndim
840 my_index(i) = iwork2(iwork1(i))
841 END DO
842 ! Iterative sorting
843 DO k = 2, nvar
844 wtmp(1:k - 1) = pot_par(1, 1:k - 1)
845 istart = 1
846 at_least_one = .false.
847 DO j = 1, ndim
848 rwork(j) = pot_par(j, k)
849 IF (all(pot_par(j, 1:k - 1) == wtmp(1:k - 1))) cycle
850 iend = j - 1
851 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 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 istart = j
860 END DO
861 iend = ndim
862 idim = iend - istart + 1
863 CALL sort(rwork(istart:iend), idim, iwork1(istart:iend))
864 iwork1(istart:iend) = iwork1(istart:iend) - 1 + istart
865 IF (idim /= 1) at_least_one = .true.
866 pot_par(:, k) = rwork
867 IF (.NOT. at_least_one) EXIT
868 ! Sort other components
869 DO j = k + 1, nvar
870 rwork(:) = pot_par(:, j)
871 DO i = 1, ndim
872 pot_par(i, j) = rwork(iwork1(i))
873 END DO
874 END DO
875 iwork2(:) = my_index
876 DO i = 1, ndim
877 my_index(i) = iwork2(iwork1(i))
878 END DO
879 END DO
880 DEALLOCATE (wtmp)
881 DEALLOCATE (iwork1)
882 DEALLOCATE (iwork2)
883 DEALLOCATE (rwork)
884 !
885 ! Let's determine the number of unique potentials and tag them
886 !
887 ALLOCATE (cwork(nvar))
888 cwork(:) = pot_par(1, :)
889 locij = my_index(1)
890 CALL get_indexes(locij, ntype, tmpij0)
891 istart = 1
892 DO j = 1, ndim
893 ! Special cases for EAM and IPBV
894 locij = my_index(j)
895 CALL get_indexes(locij, ntype, tmpij)
896 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 check)
903 CASE (gp_type)
904 check = .true.
905 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 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 IF (any(potparm%pot(tmpij(1), tmpij(2))%pot%set(1)%gp%parameters /= &
910 potparm%pot(tmpij0(1), tmpij0(2))%pot%set(1)%gp%parameters)) check = .false.
911 END IF
912 END IF
913 IF (ASSOCIATED(potparm%pot(tmpij(1), tmpij(2))%pot%set(1)%gp%values) .AND. &
914 ASSOCIATED(potparm%pot(tmpij0(1), tmpij0(2))%pot%set(1)%gp%values)) THEN
915 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 IF (any(potparm%pot(tmpij(1), tmpij(2))%pot%set(1)%gp%values /= &
918 potparm%pot(tmpij0(1), tmpij0(2))%pot%set(1)%gp%values)) check = .false.
919 END IF
920 END IF
921 CASE default
922 check = .true.
923 END SELECT
924 IF (all(cwork == pot_par(j, :)) .AND. check) cycle
925 cwork(:) = pot_par(j, :)
926 nunique = nunique + 1
927 iend = j - 1
928 CALL set_potparm_index(potparm, my_index(istart:iend), pot_target, &
929 ntype, tmpij, atomic_kind_set, shift_cutoff, do_zbl)
930 !
931 DO i = istart, iend
932 locij = my_index(i)
933 CALL get_indexes(locij, ntype, tmpij)
934 tmp_index(tmpij(1), tmpij(2)) = nunique
935 tmp_index(tmpij(2), tmpij(1)) = nunique
936 END DO
937 istart = j
938 locij = my_index(j)
939 CALL get_indexes(locij, ntype, tmpij0)
940 END DO
941 nunique = nunique + 1
942 iend = ndim
943 CALL set_potparm_index(potparm, my_index(istart:iend), pot_target, &
944 ntype, tmpij, atomic_kind_set, shift_cutoff, do_zbl)
945 DO i = istart, iend
946 locij = my_index(i)
947 CALL get_indexes(locij, ntype, tmpij)
948 tmp_index(tmpij(1), tmpij(2)) = nunique
949 tmp_index(tmpij(2), tmpij(1)) = nunique
950 END DO
951 DEALLOCATE (cwork)
952 DEALLOCATE (pot_par)
953 ELSE
954 nunique = nunique + 1
955 CALL set_potparm_index(potparm, my_index, pot_target, ntype, tmpij, &
956 atomic_kind_set, shift_cutoff, do_zbl)
957 END IF
958 DEALLOCATE (my_index)
959 END DO
960 ! Multiple defined potential
961 n = 0
962 DO i = 1, ntype
963 DO j = 1, i
964 n = n + 1
965 IF (SIZE(potparm%pot(i, j)%pot%type) == 1) cycle
966 nunique = nunique + 1
967 tmp_index(i, j) = nunique
968 tmp_index(j, i) = nunique
969 !
970 CALL set_potparm_index(potparm, (/n/), multi_type, ntype, tmpij, &
971 atomic_kind_set, shift_cutoff, do_zbl)
972 END DO
973 END DO
974 ! Concluding the postprocess..
975 ALLOCATE (spline_env)
976 CALL spline_env_create(spline_env, ntype, nunique)
977 spline_env%spltab = tmp_index
978 DEALLOCATE (tmp_index)
979 CALL timestop(handle)
980 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 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 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 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 CALL timeset(routinen, handle)
1017
1018 NULLIFY (pot, pot_ref)
1019 nvalues = SIZE(my_index)
1020 IF ((pot_target == lj_type) .OR. (pot_target == lj_charmm_type)) THEN
1021 ALLOCATE (sigma6(nvalues))
1022 ALLOCATE (wrk(nvalues))
1023 min_sigma6 = huge(0.0_dp)
1024 m_epsilon = -huge(0.0_dp)
1025 DO i = 1, nvalues
1026 value = my_index(i)
1027 CALL get_indexes(value, ntype, tmpij)
1028 pot => potparm%pot(tmpij(1), tmpij(2))%pot
1029 ! Preliminary check..
1030 check = SIZE(pot%type) == 1
1031 cpassert(check)
1032
1033 sigma6(i) = pot%set(1)%lj%sigma6
1034 l_epsilon = pot%set(1)%lj%epsilon
1035 IF (sigma6(i) /= 0.0_dp) min_sigma6 = min(min_sigma6, sigma6(i))
1036 IF (sigma6(i) == 0.0_dp) sigma6(i) = -huge(0.0_dp)
1037 IF (l_epsilon /= 0.0_dp) m_epsilon = max(m_epsilon, l_epsilon)
1038 END DO
1039 CALL sort(sigma6, nvalues, wrk)
1040 min_val = my_index(wrk(nvalues))
1041 m_sigma6 = sigma6(nvalues)
1042 ! In case there are only zeros.. let's consider them properly..
1043 IF (m_sigma6 == -huge(0.0_dp)) m_sigma6 = 1.0_dp
1044 IF (m_epsilon == -huge(0.0_dp)) m_epsilon = 0.0_dp
1045 IF (min_sigma6 == huge(0.0_dp)) min_sigma6 = 0.0_dp
1046 DEALLOCATE (sigma6)
1047 DEALLOCATE (wrk)
1048 ELSE
1049 min_val = minval(my_index(:))
1050 END IF
1051 CALL get_indexes(min_val, ntype, tmpij)
1052 tmpij_out = tmpij
1053 pot => potparm%pot(tmpij(1), tmpij(2))%pot
1054 pot%undef = .true.
1055 IF (shift_cutoff) THEN
1056 hicut0 = sqrt(pot%rcutsq)
1057 IF (abs(hicut0) <= min_hicut_value) hicut0 = default_hicut_value
1058 END IF
1059 CALL init_genpot(potparm, ntype)
1060
1061 DO i = 1, nvalues
1062 value = my_index(i)
1063 CALL get_indexes(value, ntype, tmpij)
1064 pot => potparm%pot(tmpij(1), tmpij(2))%pot
1065 CALL spline_factor_create(pot%spl_f)
1066 pot%spl_f%rcutsq_f = 1.0_dp
1067 pot%spl_f%rscale = 1.0_dp
1068 pot%spl_f%fscale = 1.0_dp
1069 END DO
1070
1071 IF (any(potential_single_allocation == pot_target)) THEN
1072 DO i = 1, nvalues
1073 value = my_index(i)
1074 CALL get_indexes(value, ntype, tmpij)
1075 pot => potparm%pot(tmpij(1), tmpij(2))%pot
1076
1077 check = SIZE(pot%type) == 1
1078 cpassert(check)
1079 ! Undef potential.. this will be used to compute the splines..
1080 IF ((pot_target == lj_type) .OR. (pot_target == lj_charmm_type)) THEN
1081 l_sigma6 = pot%set(1)%lj%sigma6
1082 l_epsilon = pot%set(1)%lj%epsilon
1083 ! Undef potential.. this will be used to compute the splines..
1084 IF (pot%undef) THEN
1085 pot%set(1)%lj%sigma6 = m_sigma6
1086 pot%set(1)%lj%sigma12 = m_sigma6**2
1087 pot%set(1)%lj%epsilon = m_epsilon
1088 END IF
1089 pot%spl_f%rscale(1) = 1.0_dp
1090 pot%spl_f%fscale(1) = 0.0_dp
1091 IF (l_sigma6*l_epsilon /= 0.0_dp) THEN
1092 pot%spl_f%rcutsq_f = (min_sigma6/m_sigma6)**(1.0_dp/3.0_dp)
1093 pot%spl_f%rscale(1) = (l_sigma6/m_sigma6)**(1.0_dp/3.0_dp)
1094 pot%spl_f%fscale(1) = l_epsilon/m_epsilon
1095 END IF
1096 END IF
1097 END DO
1098 END IF
1099
1100 DO i = 1, nvalues
1101 value = my_index(i)
1102 CALL get_indexes(value, ntype, tmpij)
1103 pot => potparm%pot(tmpij(1), tmpij(2))%pot
1104
1105 IF (do_zbl) THEN
1106 atomic_kind => atomic_kind_set(tmpij(1))
1107 CALL get_atomic_kind(atomic_kind, rcov=rcovi, z=zi)
1108 atomic_kind => atomic_kind_set(tmpij(2))
1109 CALL get_atomic_kind(atomic_kind, rcov=rcovj, z=zj)
1110 CALL zbl_matching_polinomial(pot, rcovi, rcovj, real(zi, kind=dp), &
1111 REAL(zj, kind=dp))
1112 END IF
1113 ! Derivative factors
1114 pot%spl_f%dscale = pot%spl_f%fscale/pot%spl_f%rscale
1115 ! Cutoff for the potentials on splines
1116 IF (shift_cutoff) THEN
1117 ! Cutoff NonBonded
1118 pot%spl_f%cutoff = ener_pot(pot, hicut0, 0.0_dp)
1119 END IF
1120 END DO
1121
1122 ! Handle the cutoff
1123 IF (shift_cutoff) THEN
1124 pot_ref => potparm%pot(tmpij_out(1), tmpij_out(2))%pot
1125 DO i = 1, nvalues
1126 value = my_index(i)
1127 CALL get_indexes(value, ntype, tmpij)
1128 pot => potparm%pot(tmpij(1), tmpij(2))%pot
1129 IF (value == min_val) cycle
1130 ! Cutoff NonBonded
1131 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 CALL finalizef()
1135
1136 CALL timestop(handle)
1137
1138 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 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 tmp = 0
1154 ij = huge(0)
1155 DO i = 1, ndim
1156 tmp = tmp + i
1157 IF (tmp >= inind) THEN
1158 ij(1) = i
1159 ij(2) = inind - tmp + i
1160 EXIT
1161 END IF
1162 END DO
1163 END SUBROUTINE get_indexes
1164
1165END MODULE pair_potential
1166
cp_log_handling::cp_to_string
Definition cp_log_handling.F:90
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::parsef
subroutine, public parsef(i, funcstr, var)
Parse ith function string FuncStr and compile it into bytecode.
Definition fparser.F:148
fparser::finalizef
subroutine, public finalizef()
...
Definition fparser.F:101
fparser::initf
subroutine, public initf(n)
...
Definition fparser.F:130
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::lj_charmm_type
integer, parameter, public lj_charmm_type
Definition pair_potential_types.F:36
pair_potential_types::allegro_type
integer, parameter, public allegro_type
Definition pair_potential_types.F:36
pair_potential_types::bm_type
integer, parameter, public bm_type
Definition pair_potential_types.F:36
pair_potential_types::list_pot
integer, dimension(22), parameter, public list_pot
Definition pair_potential_types.F:60
pair_potential_types::gal_type
integer, parameter, public gal_type
Definition pair_potential_types.F:36
pair_potential_types::nequip_type
integer, parameter, public nequip_type
Definition pair_potential_types.F:36
pair_potential_types::wl_type
integer, parameter, public wl_type
Definition pair_potential_types.F:36
pair_potential_types::ft_type
integer, parameter, public ft_type
Definition pair_potential_types.F:36
pair_potential_types::tab_type
integer, parameter, public tab_type
Definition pair_potential_types.F:36
pair_potential_types::ftd_type
integer, parameter, public ftd_type
Definition pair_potential_types.F:36
pair_potential_types::ip_type
integer, parameter, public ip_type
Definition pair_potential_types.F:36
pair_potential_types::lj_type
integer, parameter, public lj_type
Definition pair_potential_types.F:36
pair_potential_types::deepmd_type
integer, parameter, public deepmd_type
Definition pair_potential_types.F:36
pair_potential_types::nn_type
integer, parameter, public nn_type
Definition pair_potential_types.F:36
pair_potential_types::multi_type
integer, parameter, public multi_type
Definition pair_potential_types.F:36
pair_potential_types::quip_type
integer, parameter, public quip_type
Definition pair_potential_types.F:36
pair_potential_types::gp_type
integer, parameter, public gp_type
Definition pair_potential_types.F:36
pair_potential_types::siepmann_type
integer, parameter, public siepmann_type
Definition pair_potential_types.F:36
pair_potential_types::ace_type
integer, parameter, public ace_type
Definition pair_potential_types.F:36
pair_potential_types::compare_pot
subroutine, public compare_pot(pot1, pot2, compare)
compare two different potentials
Definition pair_potential_types.F:442
pair_potential_types::gw_type
integer, parameter, public gw_type
Definition pair_potential_types.F:36
pair_potential_types::b4_type
integer, parameter, public b4_type
Definition pair_potential_types.F:36
pair_potential_types::gal21_type
integer, parameter, public gal21_type
Definition pair_potential_types.F:36
pair_potential_types::potential_single_allocation
integer, dimension(2), public potential_single_allocation
Definition pair_potential_types.F:94
pair_potential_types::ea_type
integer, parameter, public ea_type
Definition pair_potential_types.F:36
pair_potential_types::tersoff_type
integer, parameter, public tersoff_type
Definition pair_potential_types.F:36
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:138
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:556
pair_potential::init_genpot
subroutine, public init_genpot(potparm, ntype)
Initialize genpot.
Definition pair_potential.F:79
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
atomic_kind_types::atomic_kind_type
Provides all information about an atomic kind.
Definition atomic_kind_types.F:45
cp_log_handling::cp_logger_type
type of a logger, at the moment it contains just a print level starting at which level it should be l...
Definition cp_log_handling.F:140
pair_potential_types::pair_potential_pp_type
Definition pair_potential_types.F:428
pair_potential_types::pair_potential_single_type
Definition pair_potential_types.F:395
splines_types::spline_data_p_type
Definition splines_types.F:55
splines_types::spline_data_type
Data-structure that holds all needed information about a specific spline interpolation.
Definition splines_types.F:40
splines_types::spline_environment_type
Definition splines_types.F:65
splines_types::spline_factor_type
Definition splines_types.F:71