Line data Source code
1 : !--------------------------------------------------------------------------------------------------!
2 : ! CP2K: A general program to perform molecular dynamics simulations !
3 : ! Copyright 2000-2026 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, siepmann_type, tab_type, &
36 : 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 36845 : 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 36845 : CALL timeset(routineN, handle)
88 :
89 36845 : NULLIFY (pot)
90 36845 : ngp = 0
91 : ! Prescreen for general potential type
92 896172 : DO i = 1, ntype ! i: first atom type
93 63539046 : DO j = 1, i ! j: second atom type
94 62642874 : pot => potparm%pot(i, j)%pot
95 126145099 : ngp = ngp + COUNT(pot%type == gp_type)
96 : END DO
97 : END DO
98 36845 : CALL initf(ngp)
99 36845 : ngp = 0
100 896172 : DO i = 1, ntype ! i: first atom type
101 63539046 : DO j = 1, i ! j: second atom type
102 62642874 : pot => potparm%pot(i, j)%pot
103 126145099 : DO k = 1, SIZE(pot%type)
104 125285772 : IF (pot%type(k) == gp_type) THEN
105 21984 : ngp = ngp + 1
106 21984 : pot%set(k)%gp%myid = ngp
107 21984 : 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 36845 : CALL timestop(handle)
113 :
114 36845 : 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 5250 : SUBROUTINE spline_nonbond_control(spline_env, potparm, atomic_kind_set, eps_spline, &
136 : max_energy, rlow_nb, emax_spline, npoints, iw, iw2, iw3, &
137 : do_zbl, 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 5250 : CALL timeset(routineN, handle)
157 :
158 5250 : n = 0
159 5250 : ncount = 0
160 5250 : ntype = SIZE(atomic_kind_set)
161 :
162 5250 : IF (iw3 > 0) THEN
163 : WRITE (iw3, "(/,T2,A,I0,A,I0,A)") &
164 2588 : "SPLINE_INFO| Generating ", (ntype*(ntype + 1))/2, " splines for "// &
165 5176 : TRIM(ADJUSTL(nonbonded_type))//" interactions "
166 : WRITE (iw3, "(T2,A,I0,A)") &
167 2588 : " Due to ", ntype, " different atomic kinds"
168 : END IF
169 5250 : CALL init_genpot(potparm, ntype)
170 : ! Real computation of splines
171 5250 : ip = 0
172 27656 : DO i = 1, ntype
173 543080 : DO j = 1, i
174 515424 : pot => potparm%pot(i, j)%pot
175 515424 : IF (iw3 > 0 .AND. iw <= 0) THEN
176 248588 : IF (MOD(i*(i - 1)/2 + j, MAX(1, (ntype*(ntype + 1))/(2*10))) == 0) THEN
177 11106 : WRITE (UNIT=iw3, ADVANCE="NO", FMT='(2X,A3,I0)') '...', i*(i - 1)/2 + j
178 11106 : ip = ip + 1
179 11106 : IF (ip >= 11) THEN
180 96 : WRITE (iw3, *)
181 96 : ip = 0
182 : END IF
183 : END IF
184 : END IF
185 : ! Setup of Exclusion Types
186 515424 : pot%no_pp = .TRUE.
187 515424 : pot%no_mb = .TRUE.
188 1030856 : DO k = 1, SIZE(pot%type)
189 1001141 : SELECT CASE (pot%type(k))
190 : CASE (lj_type, lj_charmm_type, wl_type, gw_type, ft_type, ftd_type, ip_type, &
191 : b4_type, bm_type, gp_type, ea_type, allegro_type, nequip_type, tab_type, deepmd_type, ace_type)
192 485709 : pot%no_pp = .FALSE.
193 : CASE (tersoff_type)
194 118 : pot%no_mb = .FALSE.
195 : CASE (siepmann_type)
196 5 : pot%no_mb = .FALSE.
197 : CASE (gal_type)
198 1 : pot%no_mb = .FALSE.
199 : CASE (gal21_type)
200 1 : pot%no_mb = .FALSE.
201 : CASE (nn_type)
202 : ! Do nothing..
203 : CASE DEFAULT
204 : ! Never reach this point
205 515432 : CPABORT("")
206 : END SELECT
207 : ! Special case for EAM
208 515424 : SELECT CASE (pot%type(k))
209 : CASE (ea_type, nequip_type, allegro_type, deepmd_type, ace_type)
210 515432 : pot%no_mb = .FALSE.
211 : END SELECT
212 : END DO
213 :
214 : ! Starting SetUp of splines
215 515424 : IF (.NOT. pot%undef) CYCLE
216 31561 : ncount = ncount + 1
217 31561 : n = spline_env%spltab(i, j)
218 31561 : locut = rlow_nb
219 31561 : hicut0 = SQRT(pot%rcutsq)
220 31561 : IF (ABS(hicut0) <= MIN_HICUT_VALUE) hicut0 = DEFAULT_HICUT_VALUE
221 31561 : hicut = hicut0/SQRT(pot%spl_f%rcutsq_f)
222 :
223 31561 : energy_cutoff = pot%spl_f%cutoff
224 :
225 : ! Find the real locut according emax_spline
226 : CALL get_spline_cutoff(hicut, locut, found_locut, pot, do_zbl, &
227 31561 : energy_cutoff, emax_spline)
228 31561 : locut = MAX(locut*SQRT(pot%spl_f%rcutsq_f), rlow_nb)
229 :
230 : ! Real Generation of the Spline
231 31561 : npoints_spline = npoints
232 : CALL generate_spline_low(spline_env%spl_pp(n)%spl_p, npoints_spline, locut, &
233 : hicut, eps_spline, iw, iw2, i, j, n, ncount, max_energy, pot, &
234 : energy_cutoff, found_locut, do_zbl, atomic_kind_set, &
235 31561 : nonbonded_type)
236 :
237 31561 : pot%undef = .FALSE.
238 : ! Unique Spline working only for a pure LJ potential..
239 31561 : IF (SIZE(pot%type) == 1) THEN
240 94655 : IF (ANY(potential_single_allocation == pot%type(1))) THEN
241 : ! Restoring the proper values for the generating spline pot
242 4 : IF ((pot%type(1) == lj_type) .OR. (pot%type(1) == lj_charmm_type)) THEN
243 4 : pot%set(1)%lj%sigma6 = pot%set(1)%lj%sigma6*pot%spl_f%rscale(1)**3
244 4 : pot%set(1)%lj%sigma12 = pot%set(1)%lj%sigma6**2
245 4 : pot%set(1)%lj%epsilon = pot%set(1)%lj%epsilon*pot%spl_f%fscale(1)
246 : END IF
247 : END IF
248 : END IF
249 : ! Correct Cutoff...
250 85528 : IF (shift_cutoff) THEN
251 : pot%spl_f%cutoff = pot%spl_f%cutoff*pot%spl_f%fscale(1) - &
252 28041 : ener_pot(pot, hicut0, 0.0_dp)
253 : END IF
254 : END DO
255 : END DO
256 5250 : CALL finalizef()
257 :
258 5250 : IF (iw > 0) THEN
259 : WRITE (UNIT=iw, FMT='(/,T2,A,I0)') &
260 26240 : "SPLINE_INFO| Number of pair potential splines allocated: ", MAXVAL(spline_env%spltab)
261 : END IF
262 5250 : IF (iw3 > 0) THEN
263 : WRITE (UNIT=iw3, FMT='(/,T2,A,I0)') &
264 525444 : "SPLINE_INFO| Number of unique splines computed: ", MAXVAL(spline_env%spltab)
265 : END IF
266 :
267 5250 : CALL timestop(handle)
268 :
269 5250 : 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 31561 : 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 31561 : dx2 = (hicut - locut)
300 31561 : x = hicut
301 31561 : locut_found = locut
302 31561 : found_locut = .FALSE.
303 94683 : DO ilevel = 1, 2
304 63122 : dx2 = dx2/100.0_dp
305 5187834 : DO jx = 1, 100
306 5161898 : e = ener_pot(pot, x, energy_cutoff)
307 5161898 : IF (do_zbl) THEN
308 5098 : e = e + ener_zbl(pot, x)
309 : END IF
310 5161898 : IF (ABS(e) > emax_spline) THEN
311 37186 : locut_found = x
312 37186 : found_locut = .TRUE.
313 37186 : EXIT
314 : END IF
315 5150648 : x = x - dx2
316 : END DO
317 94683 : x = x + dx2
318 : END DO
319 31561 : locut = locut_found
320 :
321 31561 : 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 31561 : 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 31561 : NULLIFY (logger, spl_f)
374 63122 : logger => cp_get_default_logger()
375 :
376 31561 : CALL spline_factor_create(spl_f)
377 31561 : mfac = 5
378 31561 : IF (npoints > 0) THEN
379 : fixed_spline_points = .TRUE.
380 : ELSE
381 31557 : fixed_spline_points = .FALSE.
382 31557 : npoints = 20
383 31557 : IF (.NOT. found_locut) npoints = 2
384 : END IF
385 31561 : spline_data => spl_p(1)%spline_data
386 302539 : DO WHILE (.TRUE.)
387 334100 : CALL init_splinexy(spline_data, npoints + 1)
388 334100 : dx2 = (1.0_dp/locut**2 - 1.0_dp/hicut**2)/REAL(npoints, KIND=dp)
389 334100 : x2 = 1.0_dp/hicut**2
390 334100 : spline_data%x1 = x2
391 129608508 : DO jx = 1, npoints + 1
392 : ! jx: loop over 1/distance**2
393 129274408 : x = SQRT(1.0_dp/x2)
394 129274408 : e = ener_pot(pot, x, energy_cutoff)
395 129274408 : IF (do_zbl) THEN
396 6706340 : e = e + ener_zbl(pot, x)
397 : END IF
398 129274408 : spline_data%y(jx) = e
399 129608508 : x2 = x2 + dx2
400 : END DO
401 334100 : CALL init_spline(spline_data, dx=dx2)
402 : ! This is the check for required accuracy on spline setup
403 334100 : dx2 = (hicut - locut)/REAL(mfac*npoints + 1, KIND=dp)
404 334100 : x2 = locut + dx2
405 334100 : diffmax = -1.0_dp
406 334100 : xsav = hicut
407 : ! if a fixed number of points is requested, no check on its error
408 334100 : IF (fixed_spline_points) EXIT
409 644883538 : DO jx = 1, mfac*npoints
410 644698540 : x = x2
411 644698540 : e = ener_pot(pot, x, energy_cutoff)
412 644698540 : IF (do_zbl) THEN
413 33525290 : e = e + ener_zbl(pot, x)
414 : END IF
415 644698540 : IF (ABS(e) < max_energy) THEN
416 538831133 : xdum1 = ABS(e - potential_s(spl_p, x*x, xdum, spl_f, logger))
417 538831133 : diffmax = MAX(diffmax, xdum1)
418 538831133 : xsav = MIN(x, xsav)
419 : END IF
420 644698540 : x2 = x2 + dx2
421 644883538 : IF (x2 > hicut) EXIT
422 : END DO
423 334096 : 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 334100 : 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 31561 : 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 31561 : 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 31561 : CALL spline_factor_release(spl_f)
542 :
543 31561 : 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 5250 : 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 5250 : INTEGER, ALLOCATABLE, DIMENSION(:) :: Iwork1, Iwork2, my_index
568 5250 : INTEGER, ALLOCATABLE, DIMENSION(:, :) :: tmp_index
569 : LOGICAL :: at_least_one, check
570 5250 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: Cwork, Rwork, wtmp
571 5250 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: pot_par
572 :
573 5250 : CALL timeset(routineN, handle)
574 :
575 5250 : ntype = SIZE(atomic_kind_set)
576 27656 : DO i = 1, ntype
577 543080 : DO j = 1, i
578 537830 : potparm%pot(i, j)%pot%undef = .FALSE.
579 : END DO
580 : END DO
581 21000 : ALLOCATE (tmp_index(ntype, ntype))
582 : !
583 5250 : nunique = 0
584 1036098 : tmp_index = HUGE(0)
585 120750 : DO pot_target = MINVAL(list_pot), MAXVAL(list_pot)
586 115500 : ndim = 0
587 608432 : DO i = 1, ntype
588 11947760 : DO j = 1, i
589 11339328 : IF (SIZE(potparm%pot(i, j)%pot%type) /= 1) CYCLE
590 11832084 : IF (potparm%pot(i, j)%pot%type(1) == pot_target) THEN
591 515416 : tmp_index(i, j) = 1
592 515416 : tmp_index(j, i) = 1
593 515416 : ndim = ndim + 1
594 : END IF
595 : END DO
596 : END DO
597 115500 : IF (ndim == 0) CYCLE ! No potential of this kind found
598 6155 : 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 (nequip_type)
609 2 : nvar = 1 + nvar
610 : CASE (allegro_type)
611 2 : nvar = 1 + nvar
612 : CASE (ace_type)
613 6 : nvar = 2 + nvar
614 : CASE (deepmd_type)
615 2 : nvar = 2 + nvar
616 : CASE (ft_type)
617 4 : nvar = 4 + nvar
618 : CASE (ftd_type)
619 18 : nvar = 6 + nvar
620 : CASE (ip_type)
621 260 : nvar = 3 + nvar
622 : CASE (b4_type)
623 260 : nvar = 6 + nvar
624 : CASE (bm_type)
625 6 : nvar = 9 + nvar
626 : CASE (gp_type)
627 570 : nvar = 2 + nvar
628 : CASE (tersoff_type)
629 38 : nvar = 13 + nvar
630 : CASE (siepmann_type)
631 5 : nvar = 5 + nvar
632 : CASE (gal_type)
633 1 : nvar = 12 + nvar
634 : CASE (gal21_type)
635 1 : nvar = 30 + nvar
636 : CASE (nn_type)
637 2061 : nvar = nvar
638 : CASE (tab_type)
639 8 : nvar = 4 + nvar
640 : CASE DEFAULT
641 6155 : CPABORT("")
642 : END SELECT
643 : ! Setup a table of the indexes..
644 18465 : ALLOCATE (my_index(ndim))
645 6155 : n = 0
646 6155 : nk = 0
647 35752 : DO i = 1, ntype
648 971838 : DO j = 1, i
649 936086 : n = n + 1
650 936086 : IF (SIZE(potparm%pot(i, j)%pot%type) /= 1) CYCLE
651 965679 : IF (potparm%pot(i, j)%pot%type(1) == pot_target) THEN
652 515416 : nk = nk + 1
653 515416 : my_index(nk) = n
654 : END IF
655 : END DO
656 : END DO
657 6155 : IF (nvar /= 0) THEN
658 16376 : ALLOCATE (pot_par(ndim, nvar))
659 4094 : n = 0
660 4094 : nk = 0
661 23652 : DO i = 1, ntype
662 532267 : DO j = 1, i
663 508615 : n = n + 1
664 508615 : IF (SIZE(potparm%pot(i, j)%pot%type) /= 1) CYCLE
665 528169 : IF (potparm%pot(i, j)%pot%type(1) == pot_target) THEN
666 485818 : nk = nk + 1
667 485818 : my_index(nk) = n
668 480956 : SELECT CASE (pot_target)
669 : CASE (lj_type, lj_charmm_type)
670 480956 : pot_par(nk, 1) = potparm%pot(i, j)%pot%set(1)%lj%epsilon
671 480956 : pot_par(nk, 2) = potparm%pot(i, j)%pot%set(1)%lj%sigma6
672 480956 : pot_par(nk, 3) = potparm%pot(i, j)%pot%set(1)%lj%sigma12
673 : CASE (gp_type)
674 3214 : pot_par(nk, 1) = str2id(potparm%pot(i, j)%pot%set(1)%gp%potential)
675 3214 : pot_par(nk, 2) = str2id(potparm%pot(i, j)%pot%set(1)%gp%variables)
676 : CASE (wl_type)
677 1049 : pot_par(nk, 1) = potparm%pot(i, j)%pot%set(1)%willis%a
678 1049 : pot_par(nk, 2) = potparm%pot(i, j)%pot%set(1)%willis%b
679 1049 : pot_par(nk, 3) = potparm%pot(i, j)%pot%set(1)%willis%c
680 : CASE (gw_type)
681 0 : pot_par(nk, 1) = potparm%pot(i, j)%pot%set(1)%goodwin%vr0
682 0 : pot_par(nk, 2) = potparm%pot(i, j)%pot%set(1)%goodwin%m
683 0 : pot_par(nk, 3) = potparm%pot(i, j)%pot%set(1)%goodwin%mc
684 0 : pot_par(nk, 4) = potparm%pot(i, j)%pot%set(1)%goodwin%d
685 0 : pot_par(nk, 5) = potparm%pot(i, j)%pot%set(1)%goodwin%dc
686 : CASE (ea_type)
687 20 : pot_par(nk, 1) = potparm%pot(i, j)%pot%set(1)%eam%drar
688 20 : pot_par(nk, 2) = potparm%pot(i, j)%pot%set(1)%eam%drhoar
689 20 : pot_par(nk, 3) = potparm%pot(i, j)%pot%set(1)%eam%acutal
690 20 : pot_par(nk, 4) = potparm%pot(i, j)%pot%set(1)%eam%npoints
691 : CASE (nequip_type, allegro_type)
692 : pot_par(nk, 1) = str2id( &
693 12 : TRIM(potparm%pot(i, j)%pot%set(1)%nequip%pot_file_name))
694 : CASE (ace_type)
695 : pot_par(nk, 1) = str2id( &
696 18 : TRIM(potparm%pot(i, j)%pot%set(1)%ace%ace_file_name))
697 18 : pot_par(nk, 2) = potparm%pot(i, j)%pot%set(1)%ace%atom_ace_type
698 : CASE (deepmd_type)
699 : pot_par(nk, 1) = str2id( &
700 6 : TRIM(potparm%pot(i, j)%pot%set(1)%deepmd%deepmd_file_name))
701 6 : pot_par(nk, 2) = potparm%pot(i, j)%pot%set(1)%deepmd%atom_deepmd_type
702 : CASE (ft_type)
703 12 : pot_par(nk, 1) = potparm%pot(i, j)%pot%set(1)%ft%A
704 12 : pot_par(nk, 2) = potparm%pot(i, j)%pot%set(1)%ft%B
705 12 : pot_par(nk, 3) = potparm%pot(i, j)%pot%set(1)%ft%C
706 12 : pot_par(nk, 4) = potparm%pot(i, j)%pot%set(1)%ft%D
707 : CASE (ftd_type)
708 66 : pot_par(nk, 1) = potparm%pot(i, j)%pot%set(1)%ftd%A
709 66 : pot_par(nk, 2) = potparm%pot(i, j)%pot%set(1)%ftd%B
710 66 : pot_par(nk, 3) = potparm%pot(i, j)%pot%set(1)%ftd%C
711 66 : pot_par(nk, 4) = potparm%pot(i, j)%pot%set(1)%ftd%D
712 66 : pot_par(nk, 5) = potparm%pot(i, j)%pot%set(1)%ftd%BD(1)
713 66 : pot_par(nk, 6) = potparm%pot(i, j)%pot%set(1)%ftd%BD(2)
714 : CASE (ip_type)
715 48 : pot_par(nk, 1) = potparm%pot(i, j)%pot%set(1)%ipbv%rcore
716 48 : pot_par(nk, 2) = potparm%pot(i, j)%pot%set(1)%ipbv%m
717 48 : pot_par(nk, 3) = potparm%pot(i, j)%pot%set(1)%ipbv%b
718 : CASE (b4_type)
719 260 : pot_par(nk, 1) = potparm%pot(i, j)%pot%set(1)%buck4r%a
720 260 : pot_par(nk, 2) = potparm%pot(i, j)%pot%set(1)%buck4r%b
721 260 : pot_par(nk, 3) = potparm%pot(i, j)%pot%set(1)%buck4r%c
722 260 : pot_par(nk, 4) = potparm%pot(i, j)%pot%set(1)%buck4r%r1
723 260 : pot_par(nk, 5) = potparm%pot(i, j)%pot%set(1)%buck4r%r2
724 260 : pot_par(nk, 6) = potparm%pot(i, j)%pot%set(1)%buck4r%r3
725 : CASE (bm_type)
726 10 : pot_par(nk, 1) = potparm%pot(i, j)%pot%set(1)%buckmo%f0
727 10 : pot_par(nk, 2) = potparm%pot(i, j)%pot%set(1)%buckmo%a1
728 10 : pot_par(nk, 3) = potparm%pot(i, j)%pot%set(1)%buckmo%a2
729 10 : pot_par(nk, 4) = potparm%pot(i, j)%pot%set(1)%buckmo%b1
730 10 : pot_par(nk, 5) = potparm%pot(i, j)%pot%set(1)%buckmo%b2
731 10 : pot_par(nk, 6) = potparm%pot(i, j)%pot%set(1)%buckmo%c
732 10 : pot_par(nk, 7) = potparm%pot(i, j)%pot%set(1)%buckmo%d
733 10 : pot_par(nk, 8) = potparm%pot(i, j)%pot%set(1)%buckmo%r0
734 10 : pot_par(nk, 9) = potparm%pot(i, j)%pot%set(1)%buckmo%beta
735 : CASE (tersoff_type)
736 116 : pot_par(nk, 1) = potparm%pot(i, j)%pot%set(1)%tersoff%A
737 116 : pot_par(nk, 2) = potparm%pot(i, j)%pot%set(1)%tersoff%B
738 116 : pot_par(nk, 3) = potparm%pot(i, j)%pot%set(1)%tersoff%lambda1
739 116 : pot_par(nk, 4) = potparm%pot(i, j)%pot%set(1)%tersoff%lambda2
740 116 : pot_par(nk, 5) = potparm%pot(i, j)%pot%set(1)%tersoff%alpha
741 116 : pot_par(nk, 6) = potparm%pot(i, j)%pot%set(1)%tersoff%beta
742 116 : pot_par(nk, 7) = potparm%pot(i, j)%pot%set(1)%tersoff%n
743 116 : pot_par(nk, 8) = potparm%pot(i, j)%pot%set(1)%tersoff%c
744 116 : pot_par(nk, 9) = potparm%pot(i, j)%pot%set(1)%tersoff%d
745 116 : pot_par(nk, 10) = potparm%pot(i, j)%pot%set(1)%tersoff%h
746 116 : pot_par(nk, 11) = potparm%pot(i, j)%pot%set(1)%tersoff%lambda3
747 116 : pot_par(nk, 12) = potparm%pot(i, j)%pot%set(1)%tersoff%bigR
748 116 : pot_par(nk, 13) = potparm%pot(i, j)%pot%set(1)%tersoff%bigD
749 : CASE (siepmann_type)
750 5 : pot_par(nk, 1) = potparm%pot(i, j)%pot%set(1)%siepmann%B
751 5 : pot_par(nk, 2) = potparm%pot(i, j)%pot%set(1)%siepmann%D
752 5 : pot_par(nk, 3) = potparm%pot(i, j)%pot%set(1)%siepmann%E
753 5 : pot_par(nk, 4) = potparm%pot(i, j)%pot%set(1)%siepmann%F
754 5 : pot_par(nk, 5) = potparm%pot(i, j)%pot%set(1)%siepmann%beta
755 : CASE (gal_type)
756 1 : pot_par(nk, 1) = potparm%pot(i, j)%pot%set(1)%gal%epsilon
757 1 : pot_par(nk, 2) = potparm%pot(i, j)%pot%set(1)%gal%bxy
758 1 : pot_par(nk, 3) = potparm%pot(i, j)%pot%set(1)%gal%bz
759 1 : pot_par(nk, 4) = potparm%pot(i, j)%pot%set(1)%gal%r1
760 1 : pot_par(nk, 5) = potparm%pot(i, j)%pot%set(1)%gal%r2
761 1 : pot_par(nk, 6) = potparm%pot(i, j)%pot%set(1)%gal%a1
762 1 : pot_par(nk, 7) = potparm%pot(i, j)%pot%set(1)%gal%a2
763 1 : pot_par(nk, 8) = potparm%pot(i, j)%pot%set(1)%gal%a3
764 1 : pot_par(nk, 9) = potparm%pot(i, j)%pot%set(1)%gal%a4
765 1 : pot_par(nk, 10) = potparm%pot(i, j)%pot%set(1)%gal%a
766 1 : pot_par(nk, 11) = potparm%pot(i, j)%pot%set(1)%gal%b
767 1 : pot_par(nk, 12) = potparm%pot(i, j)%pot%set(1)%gal%c
768 : CASE (gal21_type)
769 1 : pot_par(nk, 1) = potparm%pot(i, j)%pot%set(1)%gal21%epsilon1
770 1 : pot_par(nk, 2) = potparm%pot(i, j)%pot%set(1)%gal21%epsilon2
771 1 : pot_par(nk, 3) = potparm%pot(i, j)%pot%set(1)%gal21%epsilon3
772 1 : pot_par(nk, 4) = potparm%pot(i, j)%pot%set(1)%gal21%bxy1
773 1 : pot_par(nk, 5) = potparm%pot(i, j)%pot%set(1)%gal21%bxy2
774 1 : pot_par(nk, 6) = potparm%pot(i, j)%pot%set(1)%gal21%bz1
775 1 : pot_par(nk, 7) = potparm%pot(i, j)%pot%set(1)%gal21%bz2
776 1 : pot_par(nk, 8) = potparm%pot(i, j)%pot%set(1)%gal21%r1
777 1 : pot_par(nk, 9) = potparm%pot(i, j)%pot%set(1)%gal21%r2
778 1 : pot_par(nk, 10) = potparm%pot(i, j)%pot%set(1)%gal21%a11
779 1 : pot_par(nk, 11) = potparm%pot(i, j)%pot%set(1)%gal21%a12
780 1 : pot_par(nk, 12) = potparm%pot(i, j)%pot%set(1)%gal21%a13
781 1 : pot_par(nk, 13) = potparm%pot(i, j)%pot%set(1)%gal21%a21
782 1 : pot_par(nk, 14) = potparm%pot(i, j)%pot%set(1)%gal21%a22
783 1 : pot_par(nk, 15) = potparm%pot(i, j)%pot%set(1)%gal21%a23
784 1 : pot_par(nk, 16) = potparm%pot(i, j)%pot%set(1)%gal21%a31
785 1 : pot_par(nk, 17) = potparm%pot(i, j)%pot%set(1)%gal21%a32
786 1 : pot_par(nk, 18) = potparm%pot(i, j)%pot%set(1)%gal21%a33
787 1 : pot_par(nk, 19) = potparm%pot(i, j)%pot%set(1)%gal21%a41
788 1 : pot_par(nk, 20) = potparm%pot(i, j)%pot%set(1)%gal21%a42
789 1 : pot_par(nk, 21) = potparm%pot(i, j)%pot%set(1)%gal21%a43
790 1 : pot_par(nk, 22) = potparm%pot(i, j)%pot%set(1)%gal21%AO1
791 1 : pot_par(nk, 23) = potparm%pot(i, j)%pot%set(1)%gal21%AO2
792 1 : pot_par(nk, 24) = potparm%pot(i, j)%pot%set(1)%gal21%BO1
793 1 : pot_par(nk, 25) = potparm%pot(i, j)%pot%set(1)%gal21%BO2
794 1 : pot_par(nk, 26) = potparm%pot(i, j)%pot%set(1)%gal21%c
795 1 : pot_par(nk, 27) = potparm%pot(i, j)%pot%set(1)%gal21%AH1
796 1 : pot_par(nk, 28) = potparm%pot(i, j)%pot%set(1)%gal21%AH2
797 1 : pot_par(nk, 29) = potparm%pot(i, j)%pot%set(1)%gal21%BH1
798 1 : pot_par(nk, 30) = potparm%pot(i, j)%pot%set(1)%gal21%BH2
799 : CASE (tab_type)
800 24 : pot_par(nk, 1) = potparm%pot(i, j)%pot%set(1)%tab%dr
801 24 : pot_par(nk, 2) = potparm%pot(i, j)%pot%set(1)%tab%rcut
802 24 : pot_par(nk, 3) = potparm%pot(i, j)%pot%set(1)%tab%npoints
803 24 : pot_par(nk, 4) = potparm%pot(i, j)%pot%set(1)%tab%index
804 : CASE (nn_type)
805 : ! no checks
806 : CASE DEFAULT
807 485818 : CPABORT("")
808 : END SELECT
809 1448070 : IF (ANY(potential_single_allocation == pot_target)) THEN
810 37536 : pot_par(nk, :) = REAL(pot_target, KIND=dp)
811 : END IF
812 : END IF
813 : END DO
814 : END DO
815 : ! Main Sorting Loop
816 12282 : ALLOCATE (Rwork(ndim))
817 8188 : ALLOCATE (Iwork1(ndim))
818 8188 : ALLOCATE (Iwork2(ndim))
819 8188 : ALLOCATE (wtmp(nvar))
820 4094 : CALL sort(pot_par(:, 1), ndim, Iwork1)
821 : ! Sort all the other components of the potential
822 13016 : DO k = 2, nvar
823 979596 : Rwork(:) = pot_par(:, k)
824 983690 : DO i = 1, ndim
825 979596 : pot_par(i, k) = Rwork(Iwork1(i))
826 : END DO
827 : END DO
828 489912 : Iwork2(:) = my_index
829 489912 : DO i = 1, ndim
830 489912 : my_index(i) = Iwork2(Iwork1(i))
831 : END DO
832 : ! Iterative sorting
833 7695 : DO k = 2, nvar
834 13153 : wtmp(1:k - 1) = pot_par(1, 1:k - 1)
835 : istart = 1
836 : at_least_one = .FALSE.
837 969928 : DO j = 1, ndim
838 964221 : Rwork(j) = pot_par(j, k)
839 2362866 : IF (ALL(pot_par(j, 1:k - 1) == wtmp(1:k - 1))) CYCLE
840 35506 : iend = j - 1
841 90797 : wtmp(1:k - 1) = pot_par(j, 1:k - 1)
842 : ! If the ordered array has no two same consecutive elements
843 : ! does not make any sense to proceed ordering the others
844 : ! related parameters..
845 35506 : idim = iend - istart + 1
846 35506 : CALL sort(Rwork(istart:iend), idim, Iwork1(istart:iend))
847 967420 : Iwork1(istart:iend) = Iwork1(istart:iend) - 1 + istart
848 35506 : IF (idim /= 1) at_least_one = .TRUE.
849 934422 : istart = j
850 : END DO
851 5707 : iend = ndim
852 5707 : idim = iend - istart + 1
853 5707 : CALL sort(Rwork(istart:iend), idim, Iwork1(istart:iend))
854 38014 : Iwork1(istart:iend) = Iwork1(istart:iend) - 1 + istart
855 5707 : IF (idim /= 1) at_least_one = .TRUE.
856 969928 : pot_par(:, k) = Rwork
857 5707 : IF (.NOT. at_least_one) EXIT
858 : ! Sort other components
859 5360 : DO j = k + 1, nvar
860 484450 : Rwork(:) = pot_par(:, j)
861 488051 : DO i = 1, ndim
862 484450 : pot_par(i, j) = Rwork(Iwork1(i))
863 : END DO
864 : END DO
865 962290 : Iwork2(:) = my_index
866 966384 : DO i = 1, ndim
867 962290 : my_index(i) = Iwork2(Iwork1(i))
868 : END DO
869 : END DO
870 4094 : DEALLOCATE (wtmp)
871 4094 : DEALLOCATE (Iwork1)
872 4094 : DEALLOCATE (Iwork2)
873 4094 : DEALLOCATE (Rwork)
874 : !
875 : ! Let's determine the number of unique potentials and tag them
876 : !
877 8188 : ALLOCATE (Cwork(nvar))
878 17110 : Cwork(:) = pot_par(1, :)
879 4094 : locij = my_index(1)
880 4094 : CALL get_indexes(locij, ntype, tmpij0)
881 4094 : istart = 1
882 489912 : DO j = 1, ndim
883 : ! Special cases for EAM and IPBV
884 485818 : locij = my_index(j)
885 485818 : CALL get_indexes(locij, ntype, tmpij)
886 68 : SELECT CASE (pot_target)
887 : CASE (ea_type, ip_type)
888 : ! check the array components
889 : CALL compare_pot(potparm%pot(tmpij(1), tmpij(2))%pot, &
890 : potparm%pot(tmpij0(1), tmpij0(2))%pot, &
891 3282 : check)
892 : CASE (gp_type)
893 3214 : check = .TRUE.
894 3214 : IF (ASSOCIATED(potparm%pot(tmpij(1), tmpij(2))%pot%set(1)%gp%parameters) .AND. &
895 : ASSOCIATED(potparm%pot(tmpij0(1), tmpij0(2))%pot%set(1)%gp%parameters)) THEN
896 3214 : IF (SIZE(potparm%pot(tmpij(1), tmpij(2))%pot%set(1)%gp%parameters) == &
897 : SIZE(potparm%pot(tmpij0(1), tmpij0(2))%pot%set(1)%gp%parameters)) THEN
898 12652 : IF (ANY(potparm%pot(tmpij(1), tmpij(2))%pot%set(1)%gp%parameters /= &
899 0 : potparm%pot(tmpij0(1), tmpij0(2))%pot%set(1)%gp%parameters)) check = .FALSE.
900 : END IF
901 : END IF
902 3214 : IF (ASSOCIATED(potparm%pot(tmpij(1), tmpij(2))%pot%set(1)%gp%values) .AND. &
903 482536 : ASSOCIATED(potparm%pot(tmpij0(1), tmpij0(2))%pot%set(1)%gp%values)) THEN
904 3214 : IF (SIZE(potparm%pot(tmpij(1), tmpij(2))%pot%set(1)%gp%values) == &
905 : SIZE(potparm%pot(tmpij0(1), tmpij0(2))%pot%set(1)%gp%values)) THEN
906 7514 : IF (ANY(potparm%pot(tmpij(1), tmpij(2))%pot%set(1)%gp%values /= &
907 2569 : potparm%pot(tmpij0(1), tmpij0(2))%pot%set(1)%gp%values)) check = .FALSE.
908 : END IF
909 : END IF
910 : CASE default
911 485818 : check = .TRUE.
912 : END SELECT
913 1880739 : IF (ALL(Cwork == pot_par(j, :)) .AND. check) CYCLE
914 99223 : Cwork(:) = pot_par(j, :)
915 25398 : nunique = nunique + 1
916 25398 : iend = j - 1
917 : CALL set_potparm_index(potparm, my_index(istart:iend), pot_target, &
918 25398 : ntype, tmpij, atomic_kind_set, shift_cutoff, do_zbl)
919 : !
920 496115 : DO i = istart, iend
921 470717 : locij = my_index(i)
922 470717 : CALL get_indexes(locij, ntype, tmpij)
923 470717 : tmp_index(tmpij(1), tmpij(2)) = nunique
924 496115 : tmp_index(tmpij(2), tmpij(1)) = nunique
925 : END DO
926 25398 : istart = j
927 25398 : locij = my_index(j)
928 489912 : CALL get_indexes(locij, ntype, tmpij0)
929 : END DO
930 4094 : nunique = nunique + 1
931 4094 : iend = ndim
932 : CALL set_potparm_index(potparm, my_index(istart:iend), pot_target, &
933 4094 : ntype, tmpij, atomic_kind_set, shift_cutoff, do_zbl)
934 19195 : DO i = istart, iend
935 15101 : locij = my_index(i)
936 15101 : CALL get_indexes(locij, ntype, tmpij)
937 15101 : tmp_index(tmpij(1), tmpij(2)) = nunique
938 19195 : tmp_index(tmpij(2), tmpij(1)) = nunique
939 : END DO
940 4094 : DEALLOCATE (Cwork)
941 4094 : DEALLOCATE (pot_par)
942 : ELSE
943 2061 : nunique = nunique + 1
944 : CALL set_potparm_index(potparm, my_index, pot_target, ntype, tmpij, &
945 2061 : atomic_kind_set, shift_cutoff, do_zbl)
946 : END IF
947 120750 : DEALLOCATE (my_index)
948 : END DO
949 : ! Multiple defined potential
950 : n = 0
951 27656 : DO i = 1, ntype
952 543080 : DO j = 1, i
953 515424 : n = n + 1
954 515424 : IF (SIZE(potparm%pot(i, j)%pot%type) == 1) CYCLE
955 8 : nunique = nunique + 1
956 8 : tmp_index(i, j) = nunique
957 8 : tmp_index(j, i) = nunique
958 : !
959 : CALL set_potparm_index(potparm, [n], multi_type, ntype, tmpij, &
960 537838 : atomic_kind_set, shift_cutoff, do_zbl)
961 : END DO
962 : END DO
963 : ! Concluding the postprocess..
964 5250 : ALLOCATE (spline_env)
965 5250 : CALL spline_env_create(spline_env, ntype, nunique)
966 1036098 : spline_env%spltab = tmp_index
967 5250 : DEALLOCATE (tmp_index)
968 5250 : CALL timestop(handle)
969 15750 : END SUBROUTINE get_nonbond_storage
970 :
971 : ! **************************************************************************************************
972 : !> \brief Trivial for non LJ potential.. gives back in the case of LJ
973 : !> the potparm with the smallest sigma..
974 : !> \param potparm ...
975 : !> \param my_index ...
976 : !> \param pot_target ...
977 : !> \param ntype ...
978 : !> \param tmpij_out ...
979 : !> \param atomic_kind_set ...
980 : !> \param shift_cutoff ...
981 : !> \param do_zbl ...
982 : !> \author Teodoro Laino [tlaino] 2007.06
983 : ! **************************************************************************************************
984 31561 : SUBROUTINE set_potparm_index(potparm, my_index, pot_target, ntype, tmpij_out, &
985 : atomic_kind_set, shift_cutoff, do_zbl)
986 :
987 : TYPE(pair_potential_pp_type), POINTER :: potparm
988 : INTEGER, INTENT(IN) :: my_index(:), pot_target, ntype
989 : INTEGER, INTENT(OUT) :: tmpij_out(2)
990 : TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
991 : LOGICAL, INTENT(IN) :: shift_cutoff, do_zbl
992 :
993 : CHARACTER(len=*), PARAMETER :: routineN = 'set_potparm_index'
994 :
995 : INTEGER :: handle, i, min_val, nvalues, tmpij(2), &
996 : value, zi, zj
997 31561 : INTEGER, ALLOCATABLE, DIMENSION(:) :: wrk
998 : LOGICAL :: check
999 : REAL(KIND=dp) :: hicut0, l_epsilon, l_sigma6, m_epsilon, &
1000 : m_sigma6, min_sigma6, rcovi, rcovj
1001 31561 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: sigma6
1002 : TYPE(atomic_kind_type), POINTER :: atomic_kind
1003 : TYPE(pair_potential_single_type), POINTER :: pot, pot_ref
1004 :
1005 31561 : CALL timeset(routineN, handle)
1006 :
1007 31561 : NULLIFY (pot, pot_ref)
1008 31561 : nvalues = SIZE(my_index)
1009 31561 : IF ((pot_target == lj_type) .OR. (pot_target == lj_charmm_type)) THEN
1010 74913 : ALLOCATE (sigma6(nvalues))
1011 74913 : ALLOCATE (wrk(nvalues))
1012 505927 : min_sigma6 = HUGE(0.0_dp)
1013 505927 : m_epsilon = -HUGE(0.0_dp)
1014 505927 : DO i = 1, nvalues
1015 480956 : value = my_index(i)
1016 480956 : CALL get_indexes(value, ntype, tmpij)
1017 480956 : pot => potparm%pot(tmpij(1), tmpij(2))%pot
1018 : ! Preliminary check..
1019 480956 : check = SIZE(pot%type) == 1
1020 480956 : CPASSERT(check)
1021 :
1022 480956 : sigma6(i) = pot%set(1)%lj%sigma6
1023 480956 : l_epsilon = pot%set(1)%lj%epsilon
1024 480956 : IF (sigma6(i) /= 0.0_dp) min_sigma6 = MIN(min_sigma6, sigma6(i))
1025 480956 : IF (sigma6(i) == 0.0_dp) sigma6(i) = -HUGE(0.0_dp)
1026 505927 : IF (l_epsilon /= 0.0_dp) m_epsilon = MAX(m_epsilon, l_epsilon)
1027 : END DO
1028 24971 : CALL sort(sigma6, nvalues, wrk)
1029 24971 : min_val = my_index(wrk(nvalues))
1030 24971 : m_sigma6 = sigma6(nvalues)
1031 : ! In case there are only zeros.. let's consider them properly..
1032 24971 : IF (m_sigma6 == -HUGE(0.0_dp)) m_sigma6 = 1.0_dp
1033 24971 : IF (m_epsilon == -HUGE(0.0_dp)) m_epsilon = 0.0_dp
1034 24971 : IF (min_sigma6 == HUGE(0.0_dp)) min_sigma6 = 0.0_dp
1035 24971 : DEALLOCATE (sigma6)
1036 24971 : DEALLOCATE (wrk)
1037 : ELSE
1038 41058 : min_val = MINVAL(my_index(:))
1039 : END IF
1040 31561 : CALL get_indexes(min_val, ntype, tmpij)
1041 31561 : tmpij_out = tmpij
1042 31561 : pot => potparm%pot(tmpij(1), tmpij(2))%pot
1043 31561 : pot%undef = .TRUE.
1044 31561 : IF (shift_cutoff) THEN
1045 28041 : hicut0 = SQRT(pot%rcutsq)
1046 28041 : IF (ABS(hicut0) <= MIN_HICUT_VALUE) hicut0 = DEFAULT_HICUT_VALUE
1047 : END IF
1048 31561 : CALL init_genpot(potparm, ntype)
1049 :
1050 546985 : DO i = 1, nvalues
1051 515424 : value = my_index(i)
1052 515424 : CALL get_indexes(value, ntype, tmpij)
1053 515424 : pot => potparm%pot(tmpij(1), tmpij(2))%pot
1054 515424 : CALL spline_factor_create(pot%spl_f)
1055 515424 : pot%spl_f%rcutsq_f = 1.0_dp
1056 1030848 : pot%spl_f%rscale = 1.0_dp
1057 1062409 : pot%spl_f%fscale = 1.0_dp
1058 : END DO
1059 :
1060 94679 : IF (ANY(potential_single_allocation == pot_target)) THEN
1061 9388 : DO i = 1, nvalues
1062 9384 : value = my_index(i)
1063 9384 : CALL get_indexes(value, ntype, tmpij)
1064 9384 : pot => potparm%pot(tmpij(1), tmpij(2))%pot
1065 :
1066 9384 : check = SIZE(pot%type) == 1
1067 9384 : CPASSERT(check)
1068 : ! Undef potential.. this will be used to compute the splines..
1069 9388 : IF ((pot_target == lj_type) .OR. (pot_target == lj_charmm_type)) THEN
1070 9384 : l_sigma6 = pot%set(1)%lj%sigma6
1071 9384 : l_epsilon = pot%set(1)%lj%epsilon
1072 : ! Undef potential.. this will be used to compute the splines..
1073 9384 : IF (pot%undef) THEN
1074 4 : pot%set(1)%lj%sigma6 = m_sigma6
1075 4 : pot%set(1)%lj%sigma12 = m_sigma6**2
1076 4 : pot%set(1)%lj%epsilon = m_epsilon
1077 : END IF
1078 9384 : pot%spl_f%rscale(1) = 1.0_dp
1079 9384 : pot%spl_f%fscale(1) = 0.0_dp
1080 9384 : IF (l_sigma6*l_epsilon /= 0.0_dp) THEN
1081 9384 : pot%spl_f%rcutsq_f = (min_sigma6/m_sigma6)**(1.0_dp/3.0_dp)
1082 9384 : pot%spl_f%rscale(1) = (l_sigma6/m_sigma6)**(1.0_dp/3.0_dp)
1083 9384 : pot%spl_f%fscale(1) = l_epsilon/m_epsilon
1084 : END IF
1085 : END IF
1086 : END DO
1087 : END IF
1088 :
1089 546985 : DO i = 1, nvalues
1090 515424 : value = my_index(i)
1091 515424 : CALL get_indexes(value, ntype, tmpij)
1092 515424 : pot => potparm%pot(tmpij(1), tmpij(2))%pot
1093 :
1094 515424 : IF (do_zbl) THEN
1095 48 : atomic_kind => atomic_kind_set(tmpij(1))
1096 48 : CALL get_atomic_kind(atomic_kind, rcov=rcovi, z=zi)
1097 48 : atomic_kind => atomic_kind_set(tmpij(2))
1098 48 : CALL get_atomic_kind(atomic_kind, rcov=rcovj, z=zj)
1099 : CALL zbl_matching_polinomial(pot, rcovi, rcovj, REAL(zi, KIND=dp), &
1100 48 : REAL(zj, KIND=dp))
1101 : END IF
1102 : ! Derivative factors
1103 1030848 : pot%spl_f%dscale = pot%spl_f%fscale/pot%spl_f%rscale
1104 : ! Cutoff for the potentials on splines
1105 546985 : IF (shift_cutoff) THEN
1106 : ! Cutoff NonBonded
1107 115004 : pot%spl_f%cutoff = ener_pot(pot, hicut0, 0.0_dp)
1108 : END IF
1109 : END DO
1110 :
1111 : ! Handle the cutoff
1112 31561 : IF (shift_cutoff) THEN
1113 28041 : pot_ref => potparm%pot(tmpij_out(1), tmpij_out(2))%pot
1114 143045 : DO i = 1, nvalues
1115 115004 : value = my_index(i)
1116 115004 : CALL get_indexes(value, ntype, tmpij)
1117 115004 : pot => potparm%pot(tmpij(1), tmpij(2))%pot
1118 115004 : IF (value == min_val) CYCLE
1119 : ! Cutoff NonBonded
1120 143045 : pot%spl_f%cutoff = pot_ref%spl_f%cutoff*pot%spl_f%fscale(1) - pot%spl_f%cutoff
1121 : END DO
1122 : END IF
1123 31561 : CALL finalizef()
1124 :
1125 31561 : CALL timestop(handle)
1126 :
1127 31561 : END SUBROUTINE set_potparm_index
1128 :
1129 : ! **************************************************************************************************
1130 : !> \brief Gives back the indices of the matrix w.r.t. the collective array index
1131 : !> \param Inind ...
1132 : !> \param ndim ...
1133 : !> \param ij ...
1134 : !> \author Teodoro Laino [tlaino] 2006.05
1135 : ! **************************************************************************************************
1136 2668881 : SUBROUTINE get_indexes(Inind, ndim, ij)
1137 : INTEGER, INTENT(IN) :: Inind, ndim
1138 : INTEGER, DIMENSION(2), INTENT(OUT) :: ij
1139 :
1140 : INTEGER :: i, tmp
1141 :
1142 2668881 : tmp = 0
1143 8006643 : ij = HUGE(0)
1144 355422173 : DO i = 1, ndim
1145 355422173 : tmp = tmp + i
1146 355422173 : IF (tmp >= Inind) THEN
1147 2668881 : ij(1) = i
1148 2668881 : ij(2) = Inind - tmp + i
1149 2668881 : EXIT
1150 : END IF
1151 : END DO
1152 2668881 : END SUBROUTINE get_indexes
1153 :
1154 : END MODULE pair_potential
1155 :
|
