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