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, mace_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 36817 : 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 36817 : CALL timeset(routineN, handle)
88 :
89 36817 : NULLIFY (pot)
90 36817 : ngp = 0
91 : ! Prescreen for general potential type
92 896064 : DO i = 1, ntype ! i: first atom type
93 63538806 : DO j = 1, i ! j: second atom type
94 62642742 : pot => potparm%pot(i, j)%pot
95 126144755 : ngp = ngp + COUNT(pot%type == gp_type)
96 : END DO
97 : END DO
98 36817 : CALL initf(ngp)
99 36817 : ngp = 0
100 896064 : DO i = 1, ntype ! i: first atom type
101 63538806 : DO j = 1, i ! j: second atom type
102 62642742 : pot => potparm%pot(i, j)%pot
103 126144755 : DO k = 1, SIZE(pot%type)
104 125285508 : IF (pot%type(k) == gp_type) THEN
105 21992 : ngp = ngp + 1
106 21992 : pot%set(k)%gp%myid = ngp
107 21992 : 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 36817 : CALL timestop(handle)
113 :
114 36817 : 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 5242 : 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 5242 : CALL timeset(routineN, handle)
157 :
158 : n = 0
159 5242 : ncount = 0
160 5242 : ntype = SIZE(atomic_kind_set)
161 :
162 5242 : 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 5242 : CALL init_genpot(potparm, ntype)
170 : ! Real computation of splines
171 5242 : ip = 0
172 27620 : DO i = 1, ntype
173 542996 : DO j = 1, i
174 515376 : pot => potparm%pot(i, j)%pot
175 515376 : IF (iw3 > 0 .AND. iw <= 0) THEN
176 248566 : IF (MOD(i*(i - 1)/2 + j, MAX(1, (ntype*(ntype + 1))/(2*10))) == 0) THEN
177 11084 : WRITE (UNIT=iw3, ADVANCE="NO", FMT='(2X,A3,I0)') '...', i*(i - 1)/2 + j
178 11084 : ip = ip + 1
179 11084 : 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 515376 : pot%no_pp = .TRUE.
187 515376 : pot%no_mb = .TRUE.
188 1030760 : DO k = 1, SIZE(pot%type)
189 1001085 : 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, mace_type, tab_type, &
192 : deepmd_type, ace_type)
193 485701 : pot%no_pp = .FALSE.
194 : CASE (tersoff_type)
195 116 : pot%no_mb = .FALSE.
196 : CASE (siepmann_type)
197 5 : pot%no_mb = .FALSE.
198 : CASE (gal_type)
199 1 : pot%no_mb = .FALSE.
200 : CASE (gal21_type)
201 1 : pot%no_mb = .FALSE.
202 : CASE (nn_type)
203 : ! Do nothing..
204 : CASE DEFAULT
205 : ! Never reach this point
206 515384 : CPABORT("Unknown potential type for spline_nonbond_control")
207 : END SELECT
208 : ! Special case for EAM
209 515376 : SELECT CASE (pot%type(k))
210 : CASE (ea_type, nequip_type, allegro_type, mace_type, deepmd_type, ace_type)
211 515384 : pot%no_mb = .FALSE.
212 : END SELECT
213 : END DO
214 :
215 : ! Starting SetUp of splines
216 515376 : IF (.NOT. pot%undef) CYCLE
217 31541 : ncount = ncount + 1
218 31541 : n = spline_env%spltab(i, j)
219 31541 : locut = rlow_nb
220 31541 : hicut0 = SQRT(pot%rcutsq)
221 31541 : IF (ABS(hicut0) <= MIN_HICUT_VALUE) hicut0 = DEFAULT_HICUT_VALUE
222 31541 : hicut = hicut0/SQRT(pot%spl_f%rcutsq_f)
223 :
224 31541 : energy_cutoff = pot%spl_f%cutoff
225 :
226 : ! Find the real locut according emax_spline
227 : CALL get_spline_cutoff(hicut, locut, found_locut, pot, do_zbl, &
228 31541 : energy_cutoff, emax_spline)
229 31541 : locut = MAX(locut*SQRT(pot%spl_f%rcutsq_f), rlow_nb)
230 :
231 : ! Real Generation of the Spline
232 31541 : npoints_spline = npoints
233 : CALL generate_spline_low(spline_env%spl_pp(n)%spl_p, npoints_spline, locut, &
234 : hicut, eps_spline, iw, iw2, i, j, n, ncount, max_energy, pot, &
235 : energy_cutoff, found_locut, do_zbl, atomic_kind_set, &
236 31541 : nonbonded_type)
237 :
238 31541 : pot%undef = .FALSE.
239 : ! Unique Spline working only for a pure LJ potential..
240 31541 : IF (SIZE(pot%type) == 1) THEN
241 94595 : IF (ANY(potential_single_allocation == pot%type(1))) THEN
242 : ! Restoring the proper values for the generating spline pot
243 4 : IF ((pot%type(1) == lj_type) .OR. (pot%type(1) == lj_charmm_type)) THEN
244 4 : pot%set(1)%lj%sigma6 = pot%set(1)%lj%sigma6*pot%spl_f%rscale(1)**3
245 4 : pot%set(1)%lj%sigma12 = pot%set(1)%lj%sigma6**2
246 4 : pot%set(1)%lj%epsilon = pot%set(1)%lj%epsilon*pot%spl_f%fscale(1)
247 : END IF
248 : END IF
249 : END IF
250 : ! Correct Cutoff...
251 85460 : IF (shift_cutoff) THEN
252 : pot%spl_f%cutoff = pot%spl_f%cutoff*pot%spl_f%fscale(1) - &
253 28021 : ener_pot(pot, hicut0, 0.0_dp)
254 : END IF
255 : END DO
256 : END DO
257 5242 : CALL finalizef()
258 :
259 5242 : IF (iw > 0) THEN
260 : WRITE (UNIT=iw, FMT='(/,T2,A,I0)') &
261 26258 : "SPLINE_INFO| Number of pair potential splines allocated: ", MAXVAL(spline_env%spltab)
262 : END IF
263 5242 : IF (iw3 > 0) THEN
264 : WRITE (UNIT=iw3, FMT='(/,T2,A,I0)') &
265 525412 : "SPLINE_INFO| Number of unique splines computed: ", MAXVAL(spline_env%spltab)
266 : END IF
267 :
268 5242 : CALL timestop(handle)
269 :
270 5242 : END SUBROUTINE spline_nonbond_control
271 :
272 : ! **************************************************************************************************
273 : !> \brief Finds the cutoff for the generation of the spline
274 : !> In a two pass approach, first with low resolution, refine in a second iteration
275 : !> \param hicut ...
276 : !> \param locut ...
277 : !> \param found_locut ...
278 : !> \param pot ...
279 : !> \param do_zbl ...
280 : !> \param energy_cutoff ...
281 : !> \param emax_spline ...
282 : !> \par History
283 : !> Splitting in order to make some season cleaning..
284 : !> \author Teodoro Laino [tlaino] 2007.06
285 : ! **************************************************************************************************
286 31541 : SUBROUTINE get_spline_cutoff(hicut, locut, found_locut, pot, do_zbl, &
287 : energy_cutoff, emax_spline)
288 :
289 : REAL(KIND=dp), INTENT(IN) :: hicut
290 : REAL(KIND=dp), INTENT(INOUT) :: locut
291 : LOGICAL, INTENT(OUT) :: found_locut
292 : TYPE(pair_potential_single_type), OPTIONAL, &
293 : POINTER :: pot
294 : LOGICAL, INTENT(IN) :: do_zbl
295 : REAL(KIND=dp), INTENT(IN) :: energy_cutoff, emax_spline
296 :
297 : INTEGER :: ilevel, jx
298 : REAL(KIND=dp) :: dx2, e, locut_found, x
299 :
300 31541 : dx2 = (hicut - locut)
301 31541 : x = hicut
302 31541 : locut_found = locut
303 31541 : found_locut = .FALSE.
304 94623 : DO ilevel = 1, 2
305 63082 : dx2 = dx2/100.0_dp
306 5182580 : DO jx = 1, 100
307 5156716 : e = ener_pot(pot, x, energy_cutoff)
308 5156716 : IF (do_zbl) THEN
309 5098 : e = e + ener_zbl(pot, x)
310 : END IF
311 5156716 : IF (ABS(e) > emax_spline) THEN
312 37218 : locut_found = x
313 37218 : found_locut = .TRUE.
314 37218 : EXIT
315 : END IF
316 5145362 : x = x - dx2
317 : END DO
318 94623 : x = x + dx2
319 : END DO
320 31541 : locut = locut_found
321 :
322 31541 : END SUBROUTINE get_spline_cutoff
323 :
324 : ! **************************************************************************************************
325 : !> \brief Real Generation of spline..
326 : !> \param spl_p ...
327 : !> \param npoints ...
328 : !> \param locut ...
329 : !> \param hicut ...
330 : !> \param eps_spline ...
331 : !> \param iw ...
332 : !> \param iw2 ...
333 : !> \param i ...
334 : !> \param j ...
335 : !> \param n ...
336 : !> \param ncount ...
337 : !> \param max_energy ...
338 : !> \param pot ...
339 : !> \param energy_cutoff ...
340 : !> \param found_locut ...
341 : !> \param do_zbl ...
342 : !> \param atomic_kind_set ...
343 : !> \param nonbonded_type ...
344 : !> \par History
345 : !> Splitting in order to make some season cleaning..
346 : !> \author Teodoro Laino [tlaino] 2007.06
347 : ! **************************************************************************************************
348 31541 : SUBROUTINE generate_spline_low(spl_p, npoints, locut, hicut, eps_spline, &
349 : iw, iw2, i, j, n, ncount, max_energy, pot, energy_cutoff, &
350 : found_locut, do_zbl, atomic_kind_set, nonbonded_type)
351 :
352 : TYPE(spline_data_p_type), DIMENSION(:), POINTER :: spl_p
353 : INTEGER, INTENT(INOUT) :: npoints
354 : REAL(KIND=dp), INTENT(IN) :: locut, hicut, eps_spline
355 : INTEGER, INTENT(IN) :: iw, iw2, i, j, n, ncount
356 : REAL(KIND=dp), INTENT(IN) :: max_energy
357 : TYPE(pair_potential_single_type), POINTER :: pot
358 : REAL(KIND=dp), INTENT(IN), OPTIONAL :: energy_cutoff
359 : LOGICAL, INTENT(IN) :: found_locut, do_zbl
360 : TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
361 : CHARACTER(LEN=*), INTENT(IN) :: nonbonded_type
362 :
363 : CHARACTER(LEN=2*default_string_length) :: message, tmp
364 : CHARACTER(LEN=default_path_length) :: file_name
365 : INTEGER :: ix, jx, mfac, nppa, nx, unit_number
366 : LOGICAL :: fixed_spline_points
367 : REAL(KIND=dp) :: df, dg, dh, diffmax, dx, dx2, e, &
368 : e_spline, f, g, h, r, rcut, x, x2, &
369 : xdum, xdum1, xsav
370 : TYPE(cp_logger_type), POINTER :: logger
371 : TYPE(spline_data_type), POINTER :: spline_data
372 : TYPE(spline_factor_type), POINTER :: spl_f
373 :
374 31541 : NULLIFY (logger, spl_f)
375 63082 : logger => cp_get_default_logger()
376 :
377 31541 : CALL spline_factor_create(spl_f)
378 31541 : mfac = 5
379 31541 : IF (npoints > 0) THEN
380 : fixed_spline_points = .TRUE.
381 : ELSE
382 31537 : fixed_spline_points = .FALSE.
383 31537 : npoints = 20
384 31537 : IF (.NOT. found_locut) npoints = 2
385 : END IF
386 31541 : spline_data => spl_p(1)%spline_data
387 301497 : DO WHILE (.TRUE.)
388 333038 : CALL init_splinexy(spline_data, npoints + 1)
389 333038 : dx2 = (1.0_dp/locut**2 - 1.0_dp/hicut**2)/REAL(npoints, KIND=dp)
390 333038 : x2 = 1.0_dp/hicut**2
391 333038 : spline_data%x1 = x2
392 126249042 : DO jx = 1, npoints + 1
393 : ! jx: loop over 1/distance**2
394 125916004 : x = SQRT(1.0_dp/x2)
395 125916004 : e = ener_pot(pot, x, energy_cutoff)
396 125916004 : IF (do_zbl) THEN
397 6706340 : e = e + ener_zbl(pot, x)
398 : END IF
399 125916004 : spline_data%y(jx) = e
400 126249042 : x2 = x2 + dx2
401 : END DO
402 333038 : CALL init_spline(spline_data, dx=dx2)
403 : ! This is the check for required accuracy on spline setup
404 333038 : dx2 = (hicut - locut)/REAL(mfac*npoints + 1, KIND=dp)
405 333038 : x2 = locut + dx2
406 333038 : diffmax = -1.0_dp
407 333038 : xsav = hicut
408 : ! if a fixed number of points is requested, no check on its error
409 333038 : IF (fixed_spline_points) EXIT
410 628096282 : DO jx = 1, mfac*npoints
411 627911830 : x = x2
412 627911830 : e = ener_pot(pot, x, energy_cutoff)
413 627911830 : IF (do_zbl) THEN
414 33525290 : e = e + ener_zbl(pot, x)
415 : END IF
416 627911830 : IF (ABS(e) < max_energy) THEN
417 525313411 : xdum1 = ABS(e - potential_s(spl_p, x*x, xdum, spl_f, logger))
418 525313411 : diffmax = MAX(diffmax, xdum1)
419 525313411 : xsav = MIN(x, xsav)
420 : END IF
421 627911830 : x2 = x2 + dx2
422 628096282 : IF (x2 > hicut) EXIT
423 : END DO
424 333034 : IF (npoints > MAX_POINTS) THEN
425 0 : WRITE (message, '(A,I8,A,G12.6,A)') "SPLINE_INFO| Number of points: ", npoints, &
426 0 : " obtained accuracy ", diffmax, ". MM SPLINE: no convergence on required"// &
427 0 : " accuracy (adjust EPS_SPLINE and rerun)"
428 0 : CALL cp_abort(__LOCATION__, TRIM(message))
429 : END IF
430 : ! accuracy is poor or we have found no points below max_energy, refine mesh
431 333038 : IF (diffmax > eps_spline .OR. diffmax < 0.0_dp) THEN
432 301497 : npoints = CEILING(1.2_dp*REAL(npoints, KIND=dp))
433 : ELSE
434 : EXIT
435 : END IF
436 : END DO
437 : ! Print spline info to STDOUT if requested
438 31541 : IF (iw > 0) THEN
439 : WRITE (UNIT=iw, &
440 : FMT="(/,A,I0,/,A,I0,/,A,I0,1X,I0,/,A,/,A,I0,2(/,A,ES13.6),2(/,A,2ES13.6))") &
441 4806 : " SPLINE_INFO| Spline number: ", ncount, &
442 4806 : " SPLINE_INFO| Unique spline number: ", n, &
443 4806 : " SPLINE_INFO| Atomic kind numbers: ", i, j, &
444 : " SPLINE_INFO| Atomic kind names: "//TRIM(ADJUSTL(atomic_kind_set(i)%name))//" "// &
445 4806 : TRIM(ADJUSTL(atomic_kind_set(j)%name)), &
446 4806 : " SPLINE_INFO| Number of spline points: ", npoints, &
447 4806 : " SPLINE_INFO| Requested accuracy [Hartree]: ", eps_spline, &
448 4806 : " SPLINE_INFO| Achieved accuracy [Hartree]: ", diffmax, &
449 4806 : " SPLINE_INFO| Spline range [bohr]: ", locut, hicut, &
450 9612 : " SPLINE_INFO| Spline range used to achieve accuracy [bohr]:", xsav, hicut
451 4806 : dx2 = (hicut - locut)/REAL(npoints + 1, KIND=dp)
452 4806 : x = locut + dx2
453 : WRITE (UNIT=iw, FMT='(A,ES17.9)') &
454 4806 : " SPLINE_INFO| Spline value at RMIN [Hartree]: ", potential_s(spl_p, x*x, xdum, spl_f, logger), &
455 4806 : " SPLINE_INFO| Spline value at RMAX [Hartree]: ", potential_s(spl_p, hicut*hicut, xdum, spl_f, logger), &
456 9612 : " SPLINE_INFO| Non-bonded energy cutoff [Hartree]: ", energy_cutoff
457 : END IF
458 : ! Print spline data on file if requested
459 31541 : IF (iw2 > 0) THEN
460 : ! Set increment to 200 points per Angstrom
461 64 : nppa = 200
462 64 : dx = bohr/REAL(nppa, KIND=dp)
463 64 : nx = NINT(hicut/dx)
464 64 : file_name = ""
465 64 : tmp = ADJUSTL(cp_to_string(n))
466 : WRITE (UNIT=file_name, FMT="(A,I0,A)") &
467 : TRIM(ADJUSTL(nonbonded_type))//"_SPLINE_"//TRIM(tmp)//"_"// &
468 : TRIM(ADJUSTL(atomic_kind_set(i)%name))//"_"// &
469 64 : TRIM(ADJUSTL(atomic_kind_set(j)%name))
470 : CALL open_file(file_name=file_name, &
471 : file_status="UNKNOWN", &
472 : file_form="FORMATTED", &
473 : file_action="WRITE", &
474 64 : unit_number=unit_number)
475 : WRITE (UNIT=unit_number, &
476 : 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)") &
477 64 : "# Spline number: ", ncount, &
478 64 : "# Unique spline number: ", n, &
479 64 : "# Atomic kind numbers: ", i, j, &
480 : "# Atomic kind names: "//TRIM(ADJUSTL(atomic_kind_set(i)%name))//" "// &
481 64 : TRIM(ADJUSTL(atomic_kind_set(j)%name)), &
482 64 : "# Number of spline points: ", npoints, &
483 64 : "# Requested accuracy [eV]: ", eps_spline*evolt, &
484 64 : "# Achieved accuracy [eV]: ", diffmax*evolt, &
485 64 : "# Spline range [Angstrom]: ", locut/bohr, hicut/bohr, &
486 64 : "# Spline range used to achieve accuracy [Angstrom]:", xsav/bohr, hicut/bohr, &
487 64 : "# Non-bonded energy cutoff [eV]: ", energy_cutoff*evolt, &
488 64 : "# Test spline using ", nppa, " points per Angstrom:", &
489 : "# Abscissa [Angstrom] Energy [eV] Splined energy [eV] Derivative [eV/Angstrom]"// &
490 128 : " |Energy error| [eV]"
491 64 : x = 0.0_dp
492 128296 : DO jx = 0, nx
493 128232 : IF (x > hicut) x = hicut
494 128232 : IF (x > locut) THEN
495 106526 : e = ener_pot(pot, x, energy_cutoff)
496 106526 : IF (do_zbl) e = e + ener_zbl(pot, x)
497 106526 : e_spline = potential_s(spl_p, x*x, xdum, spl_f, logger)
498 : WRITE (UNIT=unit_number, FMT="(5ES25.12)") &
499 106526 : x/bohr, e*evolt, e_spline*evolt, -bohr*x*xdum*evolt, ABS((e - e_spline)*evolt)
500 : END IF
501 128296 : x = x + dx
502 : END DO
503 64 : CALL close_file(unit_number=unit_number)
504 : !MK Write table.xvf for GROMACS 4.5.5
505 : WRITE (UNIT=file_name, FMT="(A,I0,A)") &
506 : "table_"// &
507 : TRIM(ADJUSTL(atomic_kind_set(i)%name))//"_"// &
508 64 : TRIM(ADJUSTL(atomic_kind_set(j)%name))//".xvg"
509 : CALL open_file(file_name=file_name, &
510 : file_status="UNKNOWN", &
511 : file_form="FORMATTED", &
512 : file_action="WRITE", &
513 64 : unit_number=unit_number)
514 : ! Recommended increment for dp is 0.0005 nm = 0.005 Angstrom
515 : ! which are 200 points/Angstrom
516 64 : rcut = 0.1_dp*hicut/bohr
517 64 : x = 0.0_dp
518 128296 : DO jx = 0, nx
519 128232 : IF (x > hicut) x = hicut
520 128232 : r = 0.1_dp*x/bohr ! Convert bohr to nm
521 128232 : IF (x <= locut) THEN
522 : WRITE (UNIT=unit_number, FMT="(7ES25.12)") &
523 151942 : r, (0.0_dp, ix=1, 6)
524 : ELSE
525 106526 : e_spline = potential_s(spl_p, x*x, xdum, spl_f, logger)
526 106526 : f = 1.0_dp/r
527 106526 : df = -1.0_dp/r**2
528 106526 : g = -1.0_dp/r**6 + 1.0_dp/rcut**6
529 106526 : dg = 6.0_dp/r**7
530 106526 : h = e_spline*kjmol
531 106526 : dh = -10.0_dp*bohr*x*xdum*kjmol
532 : WRITE (UNIT=unit_number, FMT="(7ES25.12)") &
533 106526 : r, f, -df, & ! r, f(r), -f'(r) => probably not used
534 106526 : g, -dg, & ! g(r), -g'(r) => not used, if C = 0
535 213052 : h, -dh ! h(r), -h'(r) => used, if A = 1
536 : END IF
537 128296 : x = x + dx
538 : END DO
539 64 : CALL close_file(unit_number=unit_number)
540 : END IF
541 :
542 31541 : CALL spline_factor_release(spl_f)
543 :
544 31541 : END SUBROUTINE generate_spline_low
545 :
546 : ! **************************************************************************************************
547 : !> \brief Prescreening of the effective bonds evaluations. linear scaling algorithm
548 : !> \param spline_env ...
549 : !> \param potparm ...
550 : !> \param atomic_kind_set ...
551 : !> \param do_zbl ...
552 : !> \param shift_cutoff ...
553 : !> \author Teodoro Laino [tlaino] 2006.05
554 : ! **************************************************************************************************
555 5242 : SUBROUTINE get_nonbond_storage(spline_env, potparm, atomic_kind_set, do_zbl, &
556 : shift_cutoff)
557 :
558 : TYPE(spline_environment_type), POINTER :: spline_env
559 : TYPE(pair_potential_pp_type), POINTER :: potparm
560 : TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
561 : LOGICAL, INTENT(IN) :: do_zbl, shift_cutoff
562 :
563 : CHARACTER(len=*), PARAMETER :: routineN = 'get_nonbond_storage'
564 :
565 : INTEGER :: handle, i, idim, iend, istart, j, k, &
566 : locij, n, ndim, nk, ntype, nunique, &
567 : nvar, pot_target, tmpij(2), tmpij0(2)
568 5242 : INTEGER, ALLOCATABLE, DIMENSION(:) :: Iwork1, Iwork2, my_index
569 5242 : INTEGER, ALLOCATABLE, DIMENSION(:, :) :: tmp_index
570 : LOGICAL :: at_least_one, check
571 5242 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: Cwork, Rwork, wtmp
572 5242 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: pot_par
573 :
574 5242 : CALL timeset(routineN, handle)
575 :
576 5242 : ntype = SIZE(atomic_kind_set)
577 27620 : DO i = 1, ntype
578 542996 : DO j = 1, i
579 537754 : potparm%pot(i, j)%pot%undef = .FALSE.
580 : END DO
581 : END DO
582 20968 : ALLOCATE (tmp_index(ntype, ntype))
583 : !
584 5242 : nunique = 0
585 1035994 : tmp_index = HUGE(0)
586 125808 : DO pot_target = MINVAL(list_pot), MAXVAL(list_pot)
587 120566 : ndim = 0
588 635260 : DO i = 1, ntype
589 12488908 : DO j = 1, i
590 11853648 : IF (SIZE(potparm%pot(i, j)%pot%type) /= 1) CYCLE
591 12368158 : IF (potparm%pot(i, j)%pot%type(1) == pot_target) THEN
592 515368 : tmp_index(i, j) = 1
593 515368 : tmp_index(j, i) = 1
594 515368 : ndim = ndim + 1
595 : END IF
596 : END DO
597 : END DO
598 120566 : IF (ndim == 0) CYCLE ! No potential of this kind found
599 6137 : nvar = 0
600 : SELECT CASE (pot_target)
601 : CASE (lj_type, lj_charmm_type)
602 : nvar = 3 + nvar
603 : CASE (wl_type)
604 0 : nvar = 3 + nvar
605 : CASE (gw_type)
606 0 : nvar = 5 + nvar
607 : CASE (ea_type)
608 12 : nvar = 4 + nvar
609 : CASE (nequip_type, mace_type)
610 4 : nvar = 1 + nvar
611 : CASE (allegro_type)
612 2 : nvar = 1 + nvar
613 : CASE (ace_type)
614 6 : nvar = 2 + nvar
615 : CASE (deepmd_type)
616 2 : nvar = 2 + nvar
617 : CASE (ft_type)
618 4 : nvar = 4 + nvar
619 : CASE (ftd_type)
620 18 : nvar = 6 + nvar
621 : CASE (ip_type)
622 250 : nvar = 3 + nvar
623 : CASE (b4_type)
624 250 : nvar = 6 + nvar
625 : CASE (bm_type)
626 6 : nvar = 9 + nvar
627 : CASE (gp_type)
628 574 : nvar = 2 + nvar
629 : CASE (tersoff_type)
630 36 : nvar = 13 + nvar
631 : CASE (siepmann_type)
632 5 : nvar = 5 + nvar
633 : CASE (gal_type)
634 1 : nvar = 12 + nvar
635 : CASE (gal21_type)
636 1 : nvar = 30 + nvar
637 : CASE (nn_type)
638 2051 : nvar = nvar
639 : CASE (tab_type)
640 8 : nvar = 4 + nvar
641 : CASE DEFAULT
642 6137 : CPABORT("Unknown potential target for get_nonbond_storage")
643 : END SELECT
644 : ! Setup a table of the indexes..
645 18411 : ALLOCATE (my_index(ndim))
646 6137 : n = 0
647 6137 : nk = 0
648 35686 : DO i = 1, ntype
649 971694 : DO j = 1, i
650 936008 : n = n + 1
651 936008 : IF (SIZE(potparm%pot(i, j)%pot%type) /= 1) CYCLE
652 965553 : IF (potparm%pot(i, j)%pot%type(1) == pot_target) THEN
653 515368 : nk = nk + 1
654 515368 : my_index(nk) = n
655 : END IF
656 : END DO
657 : END DO
658 6137 : IF (nvar /= 0) THEN
659 16344 : ALLOCATE (pot_par(ndim, nvar))
660 4086 : n = 0
661 4086 : nk = 0
662 23620 : DO i = 1, ntype
663 532195 : DO j = 1, i
664 508575 : n = n + 1
665 508575 : IF (SIZE(potparm%pot(i, j)%pot%type) /= 1) CYCLE
666 528105 : IF (potparm%pot(i, j)%pot%type(1) == pot_target) THEN
667 485808 : nk = nk + 1
668 485808 : my_index(nk) = n
669 480984 : SELECT CASE (pot_target)
670 : CASE (lj_type, lj_charmm_type)
671 480984 : pot_par(nk, 1) = potparm%pot(i, j)%pot%set(1)%lj%epsilon
672 480984 : pot_par(nk, 2) = potparm%pot(i, j)%pot%set(1)%lj%sigma6
673 480984 : pot_par(nk, 3) = potparm%pot(i, j)%pot%set(1)%lj%sigma12
674 : CASE (gp_type)
675 3218 : pot_par(nk, 1) = str2id(potparm%pot(i, j)%pot%set(1)%gp%potential)
676 3218 : pot_par(nk, 2) = str2id(potparm%pot(i, j)%pot%set(1)%gp%variables)
677 : CASE (wl_type)
678 1017 : pot_par(nk, 1) = potparm%pot(i, j)%pot%set(1)%willis%a
679 1017 : pot_par(nk, 2) = potparm%pot(i, j)%pot%set(1)%willis%b
680 1017 : pot_par(nk, 3) = potparm%pot(i, j)%pot%set(1)%willis%c
681 : CASE (gw_type)
682 0 : pot_par(nk, 1) = potparm%pot(i, j)%pot%set(1)%goodwin%vr0
683 0 : pot_par(nk, 2) = potparm%pot(i, j)%pot%set(1)%goodwin%m
684 0 : pot_par(nk, 3) = potparm%pot(i, j)%pot%set(1)%goodwin%mc
685 0 : pot_par(nk, 4) = potparm%pot(i, j)%pot%set(1)%goodwin%d
686 0 : pot_par(nk, 5) = potparm%pot(i, j)%pot%set(1)%goodwin%dc
687 : CASE (ea_type)
688 20 : pot_par(nk, 1) = potparm%pot(i, j)%pot%set(1)%eam%drar
689 20 : pot_par(nk, 2) = potparm%pot(i, j)%pot%set(1)%eam%drhoar
690 20 : pot_par(nk, 3) = potparm%pot(i, j)%pot%set(1)%eam%acutal
691 20 : pot_par(nk, 4) = potparm%pot(i, j)%pot%set(1)%eam%npoints
692 : CASE (nequip_type, allegro_type, mace_type)
693 : pot_par(nk, 1) = str2id( &
694 14 : TRIM(potparm%pot(i, j)%pot%set(1)%nequip%pot_file_name))
695 : CASE (ace_type)
696 : pot_par(nk, 1) = str2id( &
697 18 : TRIM(potparm%pot(i, j)%pot%set(1)%ace%ace_file_name))
698 18 : pot_par(nk, 2) = potparm%pot(i, j)%pot%set(1)%ace%atom_ace_type
699 : CASE (deepmd_type)
700 : pot_par(nk, 1) = str2id( &
701 6 : TRIM(potparm%pot(i, j)%pot%set(1)%deepmd%deepmd_file_name))
702 6 : pot_par(nk, 2) = potparm%pot(i, j)%pot%set(1)%deepmd%atom_deepmd_type
703 : CASE (ft_type)
704 12 : pot_par(nk, 1) = potparm%pot(i, j)%pot%set(1)%ft%A
705 12 : pot_par(nk, 2) = potparm%pot(i, j)%pot%set(1)%ft%B
706 12 : pot_par(nk, 3) = potparm%pot(i, j)%pot%set(1)%ft%C
707 12 : pot_par(nk, 4) = potparm%pot(i, j)%pot%set(1)%ft%D
708 : CASE (ftd_type)
709 66 : pot_par(nk, 1) = potparm%pot(i, j)%pot%set(1)%ftd%A
710 66 : pot_par(nk, 2) = potparm%pot(i, j)%pot%set(1)%ftd%B
711 66 : pot_par(nk, 3) = potparm%pot(i, j)%pot%set(1)%ftd%C
712 66 : pot_par(nk, 4) = potparm%pot(i, j)%pot%set(1)%ftd%D
713 66 : pot_par(nk, 5) = potparm%pot(i, j)%pot%set(1)%ftd%BD(1)
714 66 : pot_par(nk, 6) = potparm%pot(i, j)%pot%set(1)%ftd%BD(2)
715 : CASE (ip_type)
716 48 : pot_par(nk, 1) = potparm%pot(i, j)%pot%set(1)%ipbv%rcore
717 48 : pot_par(nk, 2) = potparm%pot(i, j)%pot%set(1)%ipbv%m
718 48 : pot_par(nk, 3) = potparm%pot(i, j)%pot%set(1)%ipbv%b
719 : CASE (b4_type)
720 250 : pot_par(nk, 1) = potparm%pot(i, j)%pot%set(1)%buck4r%a
721 250 : pot_par(nk, 2) = potparm%pot(i, j)%pot%set(1)%buck4r%b
722 250 : pot_par(nk, 3) = potparm%pot(i, j)%pot%set(1)%buck4r%c
723 250 : pot_par(nk, 4) = potparm%pot(i, j)%pot%set(1)%buck4r%r1
724 250 : pot_par(nk, 5) = potparm%pot(i, j)%pot%set(1)%buck4r%r2
725 250 : pot_par(nk, 6) = potparm%pot(i, j)%pot%set(1)%buck4r%r3
726 : CASE (bm_type)
727 10 : pot_par(nk, 1) = potparm%pot(i, j)%pot%set(1)%buckmo%f0
728 10 : pot_par(nk, 2) = potparm%pot(i, j)%pot%set(1)%buckmo%a1
729 10 : pot_par(nk, 3) = potparm%pot(i, j)%pot%set(1)%buckmo%a2
730 10 : pot_par(nk, 4) = potparm%pot(i, j)%pot%set(1)%buckmo%b1
731 10 : pot_par(nk, 5) = potparm%pot(i, j)%pot%set(1)%buckmo%b2
732 10 : pot_par(nk, 6) = potparm%pot(i, j)%pot%set(1)%buckmo%c
733 10 : pot_par(nk, 7) = potparm%pot(i, j)%pot%set(1)%buckmo%d
734 10 : pot_par(nk, 8) = potparm%pot(i, j)%pot%set(1)%buckmo%r0
735 10 : pot_par(nk, 9) = potparm%pot(i, j)%pot%set(1)%buckmo%beta
736 : CASE (tersoff_type)
737 114 : pot_par(nk, 1) = potparm%pot(i, j)%pot%set(1)%tersoff%A
738 114 : pot_par(nk, 2) = potparm%pot(i, j)%pot%set(1)%tersoff%B
739 114 : pot_par(nk, 3) = potparm%pot(i, j)%pot%set(1)%tersoff%lambda1
740 114 : pot_par(nk, 4) = potparm%pot(i, j)%pot%set(1)%tersoff%lambda2
741 114 : pot_par(nk, 5) = potparm%pot(i, j)%pot%set(1)%tersoff%alpha
742 114 : pot_par(nk, 6) = potparm%pot(i, j)%pot%set(1)%tersoff%beta
743 114 : pot_par(nk, 7) = potparm%pot(i, j)%pot%set(1)%tersoff%n
744 114 : pot_par(nk, 8) = potparm%pot(i, j)%pot%set(1)%tersoff%c
745 114 : pot_par(nk, 9) = potparm%pot(i, j)%pot%set(1)%tersoff%d
746 114 : pot_par(nk, 10) = potparm%pot(i, j)%pot%set(1)%tersoff%h
747 114 : pot_par(nk, 11) = potparm%pot(i, j)%pot%set(1)%tersoff%lambda3
748 114 : pot_par(nk, 12) = potparm%pot(i, j)%pot%set(1)%tersoff%bigR
749 114 : pot_par(nk, 13) = potparm%pot(i, j)%pot%set(1)%tersoff%bigD
750 : CASE (siepmann_type)
751 5 : pot_par(nk, 1) = potparm%pot(i, j)%pot%set(1)%siepmann%B
752 5 : pot_par(nk, 2) = potparm%pot(i, j)%pot%set(1)%siepmann%D
753 5 : pot_par(nk, 3) = potparm%pot(i, j)%pot%set(1)%siepmann%E
754 5 : pot_par(nk, 4) = potparm%pot(i, j)%pot%set(1)%siepmann%F
755 5 : pot_par(nk, 5) = potparm%pot(i, j)%pot%set(1)%siepmann%beta
756 : CASE (gal_type)
757 1 : pot_par(nk, 1) = potparm%pot(i, j)%pot%set(1)%gal%epsilon
758 1 : pot_par(nk, 2) = potparm%pot(i, j)%pot%set(1)%gal%bxy
759 1 : pot_par(nk, 3) = potparm%pot(i, j)%pot%set(1)%gal%bz
760 1 : pot_par(nk, 4) = potparm%pot(i, j)%pot%set(1)%gal%r1
761 1 : pot_par(nk, 5) = potparm%pot(i, j)%pot%set(1)%gal%r2
762 1 : pot_par(nk, 6) = potparm%pot(i, j)%pot%set(1)%gal%a1
763 1 : pot_par(nk, 7) = potparm%pot(i, j)%pot%set(1)%gal%a2
764 1 : pot_par(nk, 8) = potparm%pot(i, j)%pot%set(1)%gal%a3
765 1 : pot_par(nk, 9) = potparm%pot(i, j)%pot%set(1)%gal%a4
766 1 : pot_par(nk, 10) = potparm%pot(i, j)%pot%set(1)%gal%a
767 1 : pot_par(nk, 11) = potparm%pot(i, j)%pot%set(1)%gal%b
768 1 : pot_par(nk, 12) = potparm%pot(i, j)%pot%set(1)%gal%c
769 : CASE (gal21_type)
770 1 : pot_par(nk, 1) = potparm%pot(i, j)%pot%set(1)%gal21%epsilon1
771 1 : pot_par(nk, 2) = potparm%pot(i, j)%pot%set(1)%gal21%epsilon2
772 1 : pot_par(nk, 3) = potparm%pot(i, j)%pot%set(1)%gal21%epsilon3
773 1 : pot_par(nk, 4) = potparm%pot(i, j)%pot%set(1)%gal21%bxy1
774 1 : pot_par(nk, 5) = potparm%pot(i, j)%pot%set(1)%gal21%bxy2
775 1 : pot_par(nk, 6) = potparm%pot(i, j)%pot%set(1)%gal21%bz1
776 1 : pot_par(nk, 7) = potparm%pot(i, j)%pot%set(1)%gal21%bz2
777 1 : pot_par(nk, 8) = potparm%pot(i, j)%pot%set(1)%gal21%r1
778 1 : pot_par(nk, 9) = potparm%pot(i, j)%pot%set(1)%gal21%r2
779 1 : pot_par(nk, 10) = potparm%pot(i, j)%pot%set(1)%gal21%a11
780 1 : pot_par(nk, 11) = potparm%pot(i, j)%pot%set(1)%gal21%a12
781 1 : pot_par(nk, 12) = potparm%pot(i, j)%pot%set(1)%gal21%a13
782 1 : pot_par(nk, 13) = potparm%pot(i, j)%pot%set(1)%gal21%a21
783 1 : pot_par(nk, 14) = potparm%pot(i, j)%pot%set(1)%gal21%a22
784 1 : pot_par(nk, 15) = potparm%pot(i, j)%pot%set(1)%gal21%a23
785 1 : pot_par(nk, 16) = potparm%pot(i, j)%pot%set(1)%gal21%a31
786 1 : pot_par(nk, 17) = potparm%pot(i, j)%pot%set(1)%gal21%a32
787 1 : pot_par(nk, 18) = potparm%pot(i, j)%pot%set(1)%gal21%a33
788 1 : pot_par(nk, 19) = potparm%pot(i, j)%pot%set(1)%gal21%a41
789 1 : pot_par(nk, 20) = potparm%pot(i, j)%pot%set(1)%gal21%a42
790 1 : pot_par(nk, 21) = potparm%pot(i, j)%pot%set(1)%gal21%a43
791 1 : pot_par(nk, 22) = potparm%pot(i, j)%pot%set(1)%gal21%AO1
792 1 : pot_par(nk, 23) = potparm%pot(i, j)%pot%set(1)%gal21%AO2
793 1 : pot_par(nk, 24) = potparm%pot(i, j)%pot%set(1)%gal21%BO1
794 1 : pot_par(nk, 25) = potparm%pot(i, j)%pot%set(1)%gal21%BO2
795 1 : pot_par(nk, 26) = potparm%pot(i, j)%pot%set(1)%gal21%c
796 1 : pot_par(nk, 27) = potparm%pot(i, j)%pot%set(1)%gal21%AH1
797 1 : pot_par(nk, 28) = potparm%pot(i, j)%pot%set(1)%gal21%AH2
798 1 : pot_par(nk, 29) = potparm%pot(i, j)%pot%set(1)%gal21%BH1
799 1 : pot_par(nk, 30) = potparm%pot(i, j)%pot%set(1)%gal21%BH2
800 : CASE (tab_type)
801 24 : pot_par(nk, 1) = potparm%pot(i, j)%pot%set(1)%tab%dr
802 24 : pot_par(nk, 2) = potparm%pot(i, j)%pot%set(1)%tab%rcut
803 24 : pot_par(nk, 3) = potparm%pot(i, j)%pot%set(1)%tab%npoints
804 24 : pot_par(nk, 4) = potparm%pot(i, j)%pot%set(1)%tab%index
805 : CASE (nn_type)
806 : ! no checks
807 : CASE DEFAULT
808 485808 : CPABORT("Unknown potential target for get_nonbond_storage")
809 : END SELECT
810 1448040 : IF (ANY(potential_single_allocation == pot_target)) THEN
811 37536 : pot_par(nk, :) = REAL(pot_target, KIND=dp)
812 : END IF
813 : END IF
814 : END DO
815 : END DO
816 : ! Main Sorting Loop
817 12258 : ALLOCATE (Rwork(ndim))
818 8172 : ALLOCATE (Iwork1(ndim))
819 8172 : ALLOCATE (Iwork2(ndim))
820 8172 : ALLOCATE (wtmp(nvar))
821 4086 : CALL sort(pot_par(:, 1), ndim, Iwork1)
822 : ! Sort all the other components of the potential
823 12934 : DO k = 2, nvar
824 979444 : Rwork(:) = pot_par(:, k)
825 983530 : DO i = 1, ndim
826 979444 : pot_par(i, k) = Rwork(Iwork1(i))
827 : END DO
828 : END DO
829 489894 : Iwork2(:) = my_index
830 489894 : DO i = 1, ndim
831 489894 : my_index(i) = Iwork2(Iwork1(i))
832 : END DO
833 : ! Iterative sorting
834 7687 : DO k = 2, nvar
835 13133 : wtmp(1:k - 1) = pot_par(1, 1:k - 1)
836 : istart = 1
837 : at_least_one = .FALSE.
838 969906 : DO j = 1, ndim
839 964209 : Rwork(j) = pot_par(j, k)
840 2362844 : IF (ALL(pot_par(j, 1:k - 1) == wtmp(1:k - 1))) CYCLE
841 35504 : iend = j - 1
842 90793 : wtmp(1:k - 1) = pot_par(j, 1:k - 1)
843 : ! If the ordered array has no two same consecutive elements
844 : ! does not make any sense to proceed ordering the others
845 : ! related parameters..
846 35504 : idim = iend - istart + 1
847 35504 : CALL sort(Rwork(istart:iend), idim, Iwork1(istart:iend))
848 967416 : Iwork1(istart:iend) = Iwork1(istart:iend) - 1 + istart
849 35504 : IF (idim /= 1) at_least_one = .TRUE.
850 934402 : istart = j
851 : END DO
852 5697 : iend = ndim
853 5697 : idim = iend - istart + 1
854 5697 : CALL sort(Rwork(istart:iend), idim, Iwork1(istart:iend))
855 37994 : Iwork1(istart:iend) = Iwork1(istart:iend) - 1 + istart
856 5697 : IF (idim /= 1) at_least_one = .TRUE.
857 969906 : pot_par(:, k) = Rwork
858 5697 : IF (.NOT. at_least_one) EXIT
859 : ! Sort other components
860 5360 : DO j = k + 1, nvar
861 484450 : Rwork(:) = pot_par(:, j)
862 488051 : DO i = 1, ndim
863 484450 : pot_par(i, j) = Rwork(Iwork1(i))
864 : END DO
865 : END DO
866 962290 : Iwork2(:) = my_index
867 966376 : DO i = 1, ndim
868 962290 : my_index(i) = Iwork2(Iwork1(i))
869 : END DO
870 : END DO
871 4086 : DEALLOCATE (wtmp)
872 4086 : DEALLOCATE (Iwork1)
873 4086 : DEALLOCATE (Iwork2)
874 4086 : DEALLOCATE (Rwork)
875 : !
876 : ! Let's determine the number of unique potentials and tag them
877 : !
878 8172 : ALLOCATE (Cwork(nvar))
879 17020 : Cwork(:) = pot_par(1, :)
880 4086 : locij = my_index(1)
881 4086 : CALL get_indexes(locij, ntype, tmpij0)
882 4086 : istart = 1
883 489894 : DO j = 1, ndim
884 : ! Special cases for EAM and IPBV
885 485808 : locij = my_index(j)
886 485808 : CALL get_indexes(locij, ntype, tmpij)
887 68 : SELECT CASE (pot_target)
888 : CASE (ea_type, ip_type)
889 : ! check the array components
890 : CALL compare_pot(potparm%pot(tmpij(1), tmpij(2))%pot, &
891 : potparm%pot(tmpij0(1), tmpij0(2))%pot, &
892 3286 : check)
893 : CASE (gp_type)
894 3218 : check = .TRUE.
895 3218 : IF (ASSOCIATED(potparm%pot(tmpij(1), tmpij(2))%pot%set(1)%gp%parameters) .AND. &
896 : ASSOCIATED(potparm%pot(tmpij0(1), tmpij0(2))%pot%set(1)%gp%parameters)) THEN
897 3218 : IF (SIZE(potparm%pot(tmpij(1), tmpij(2))%pot%set(1)%gp%parameters) == &
898 : SIZE(potparm%pot(tmpij0(1), tmpij0(2))%pot%set(1)%gp%parameters)) THEN
899 12660 : IF (ANY(potparm%pot(tmpij(1), tmpij(2))%pot%set(1)%gp%parameters /= &
900 0 : potparm%pot(tmpij0(1), tmpij0(2))%pot%set(1)%gp%parameters)) check = .FALSE.
901 : END IF
902 : END IF
903 3218 : IF (ASSOCIATED(potparm%pot(tmpij(1), tmpij(2))%pot%set(1)%gp%values) .AND. &
904 482522 : ASSOCIATED(potparm%pot(tmpij0(1), tmpij0(2))%pot%set(1)%gp%values)) THEN
905 3218 : IF (SIZE(potparm%pot(tmpij(1), tmpij(2))%pot%set(1)%gp%values) == &
906 : SIZE(potparm%pot(tmpij0(1), tmpij0(2))%pot%set(1)%gp%values)) THEN
907 7522 : IF (ANY(potparm%pot(tmpij(1), tmpij(2))%pot%set(1)%gp%values /= &
908 2569 : potparm%pot(tmpij0(1), tmpij0(2))%pot%set(1)%gp%values)) check = .FALSE.
909 : END IF
910 : END IF
911 : CASE default
912 485808 : check = .TRUE.
913 : END SELECT
914 1880647 : IF (ALL(Cwork == pot_par(j, :)) .AND. check) CYCLE
915 99215 : Cwork(:) = pot_par(j, :)
916 25396 : nunique = nunique + 1
917 25396 : iend = j - 1
918 : CALL set_potparm_index(potparm, my_index(istart:iend), pot_target, &
919 25396 : ntype, tmpij, atomic_kind_set, shift_cutoff, do_zbl)
920 : !
921 496111 : DO i = istart, iend
922 470715 : locij = my_index(i)
923 470715 : CALL get_indexes(locij, ntype, tmpij)
924 470715 : tmp_index(tmpij(1), tmpij(2)) = nunique
925 496111 : tmp_index(tmpij(2), tmpij(1)) = nunique
926 : END DO
927 25396 : istart = j
928 25396 : locij = my_index(j)
929 489894 : CALL get_indexes(locij, ntype, tmpij0)
930 : END DO
931 4086 : nunique = nunique + 1
932 4086 : iend = ndim
933 : CALL set_potparm_index(potparm, my_index(istart:iend), pot_target, &
934 4086 : ntype, tmpij, atomic_kind_set, shift_cutoff, do_zbl)
935 19179 : DO i = istart, iend
936 15093 : locij = my_index(i)
937 15093 : CALL get_indexes(locij, ntype, tmpij)
938 15093 : tmp_index(tmpij(1), tmpij(2)) = nunique
939 19179 : tmp_index(tmpij(2), tmpij(1)) = nunique
940 : END DO
941 4086 : DEALLOCATE (Cwork)
942 4086 : DEALLOCATE (pot_par)
943 : ELSE
944 2051 : nunique = nunique + 1
945 : CALL set_potparm_index(potparm, my_index, pot_target, ntype, tmpij, &
946 2051 : atomic_kind_set, shift_cutoff, do_zbl)
947 : END IF
948 125808 : DEALLOCATE (my_index)
949 : END DO
950 : ! Multiple defined potential
951 : n = 0
952 27620 : DO i = 1, ntype
953 542996 : DO j = 1, i
954 515376 : n = n + 1
955 515376 : IF (SIZE(potparm%pot(i, j)%pot%type) == 1) CYCLE
956 8 : nunique = nunique + 1
957 8 : tmp_index(i, j) = nunique
958 8 : tmp_index(j, i) = nunique
959 : !
960 : CALL set_potparm_index(potparm, [n], multi_type, ntype, tmpij, &
961 537762 : atomic_kind_set, shift_cutoff, do_zbl)
962 : END DO
963 : END DO
964 : ! Concluding the postprocess..
965 5242 : ALLOCATE (spline_env)
966 5242 : CALL spline_env_create(spline_env, ntype, nunique)
967 1035994 : spline_env%spltab = tmp_index
968 5242 : DEALLOCATE (tmp_index)
969 5242 : CALL timestop(handle)
970 15726 : END SUBROUTINE get_nonbond_storage
971 :
972 : ! **************************************************************************************************
973 : !> \brief Trivial for non LJ potential.. gives back in the case of LJ
974 : !> the potparm with the smallest sigma..
975 : !> \param potparm ...
976 : !> \param my_index ...
977 : !> \param pot_target ...
978 : !> \param ntype ...
979 : !> \param tmpij_out ...
980 : !> \param atomic_kind_set ...
981 : !> \param shift_cutoff ...
982 : !> \param do_zbl ...
983 : !> \author Teodoro Laino [tlaino] 2007.06
984 : ! **************************************************************************************************
985 31541 : SUBROUTINE set_potparm_index(potparm, my_index, pot_target, ntype, tmpij_out, &
986 : atomic_kind_set, shift_cutoff, do_zbl)
987 :
988 : TYPE(pair_potential_pp_type), POINTER :: potparm
989 : INTEGER, INTENT(IN) :: my_index(:), pot_target, ntype
990 : INTEGER, INTENT(OUT) :: tmpij_out(2)
991 : TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
992 : LOGICAL, INTENT(IN) :: shift_cutoff, do_zbl
993 :
994 : CHARACTER(len=*), PARAMETER :: routineN = 'set_potparm_index'
995 :
996 : INTEGER :: handle, i, min_val, nvalues, tmpij(2), &
997 : value, zi, zj
998 31541 : INTEGER, ALLOCATABLE, DIMENSION(:) :: wrk
999 : LOGICAL :: check
1000 : REAL(KIND=dp) :: hicut0, l_epsilon, l_sigma6, m_epsilon, &
1001 : m_sigma6, min_sigma6, rcovi, rcovj
1002 31541 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: sigma6
1003 : TYPE(atomic_kind_type), POINTER :: atomic_kind
1004 : TYPE(pair_potential_single_type), POINTER :: pot, pot_ref
1005 :
1006 31541 : CALL timeset(routineN, handle)
1007 :
1008 31541 : NULLIFY (pot, pot_ref)
1009 31541 : nvalues = SIZE(my_index)
1010 31541 : IF ((pot_target == lj_type) .OR. (pot_target == lj_charmm_type)) THEN
1011 74997 : ALLOCATE (sigma6(nvalues))
1012 74997 : ALLOCATE (wrk(nvalues))
1013 505983 : min_sigma6 = HUGE(0.0_dp)
1014 505983 : m_epsilon = -HUGE(0.0_dp)
1015 505983 : DO i = 1, nvalues
1016 480984 : value = my_index(i)
1017 480984 : CALL get_indexes(value, ntype, tmpij)
1018 480984 : pot => potparm%pot(tmpij(1), tmpij(2))%pot
1019 : ! Preliminary check..
1020 480984 : check = SIZE(pot%type) == 1
1021 480984 : CPASSERT(check)
1022 :
1023 480984 : sigma6(i) = pot%set(1)%lj%sigma6
1024 480984 : l_epsilon = pot%set(1)%lj%epsilon
1025 480984 : IF (sigma6(i) /= 0.0_dp) min_sigma6 = MIN(min_sigma6, sigma6(i))
1026 480984 : IF (sigma6(i) == 0.0_dp) sigma6(i) = -HUGE(0.0_dp)
1027 505983 : IF (l_epsilon /= 0.0_dp) m_epsilon = MAX(m_epsilon, l_epsilon)
1028 : END DO
1029 24999 : CALL sort(sigma6, nvalues, wrk)
1030 24999 : min_val = my_index(wrk(nvalues))
1031 24999 : m_sigma6 = sigma6(nvalues)
1032 : ! In case there are only zeros.. let's consider them properly..
1033 24999 : IF (m_sigma6 == -HUGE(0.0_dp)) m_sigma6 = 1.0_dp
1034 24999 : IF (m_epsilon == -HUGE(0.0_dp)) m_epsilon = 0.0_dp
1035 24999 : IF (min_sigma6 == HUGE(0.0_dp)) min_sigma6 = 0.0_dp
1036 24999 : DEALLOCATE (sigma6)
1037 24999 : DEALLOCATE (wrk)
1038 : ELSE
1039 40934 : min_val = MINVAL(my_index(:))
1040 : END IF
1041 31541 : CALL get_indexes(min_val, ntype, tmpij)
1042 31541 : tmpij_out = tmpij
1043 31541 : pot => potparm%pot(tmpij(1), tmpij(2))%pot
1044 31541 : pot%undef = .TRUE.
1045 31541 : IF (shift_cutoff) THEN
1046 28021 : hicut0 = SQRT(pot%rcutsq)
1047 28021 : IF (ABS(hicut0) <= MIN_HICUT_VALUE) hicut0 = DEFAULT_HICUT_VALUE
1048 : END IF
1049 31541 : CALL init_genpot(potparm, ntype)
1050 :
1051 546917 : DO i = 1, nvalues
1052 515376 : value = my_index(i)
1053 515376 : CALL get_indexes(value, ntype, tmpij)
1054 515376 : pot => potparm%pot(tmpij(1), tmpij(2))%pot
1055 515376 : CALL spline_factor_create(pot%spl_f)
1056 515376 : pot%spl_f%rcutsq_f = 1.0_dp
1057 1030752 : pot%spl_f%rscale = 1.0_dp
1058 1062293 : pot%spl_f%fscale = 1.0_dp
1059 : END DO
1060 :
1061 94619 : IF (ANY(potential_single_allocation == pot_target)) THEN
1062 9388 : DO i = 1, nvalues
1063 9384 : value = my_index(i)
1064 9384 : CALL get_indexes(value, ntype, tmpij)
1065 9384 : pot => potparm%pot(tmpij(1), tmpij(2))%pot
1066 :
1067 9384 : check = SIZE(pot%type) == 1
1068 9384 : CPASSERT(check)
1069 : ! Undef potential.. this will be used to compute the splines..
1070 9388 : IF ((pot_target == lj_type) .OR. (pot_target == lj_charmm_type)) THEN
1071 9384 : l_sigma6 = pot%set(1)%lj%sigma6
1072 9384 : l_epsilon = pot%set(1)%lj%epsilon
1073 : ! Undef potential.. this will be used to compute the splines..
1074 9384 : IF (pot%undef) THEN
1075 4 : pot%set(1)%lj%sigma6 = m_sigma6
1076 4 : pot%set(1)%lj%sigma12 = m_sigma6**2
1077 4 : pot%set(1)%lj%epsilon = m_epsilon
1078 : END IF
1079 9384 : pot%spl_f%rscale(1) = 1.0_dp
1080 9384 : pot%spl_f%fscale(1) = 0.0_dp
1081 9384 : IF (l_sigma6*l_epsilon /= 0.0_dp) THEN
1082 9384 : pot%spl_f%rcutsq_f = (min_sigma6/m_sigma6)**(1.0_dp/3.0_dp)
1083 9384 : pot%spl_f%rscale(1) = (l_sigma6/m_sigma6)**(1.0_dp/3.0_dp)
1084 9384 : pot%spl_f%fscale(1) = l_epsilon/m_epsilon
1085 : END IF
1086 : END IF
1087 : END DO
1088 : END IF
1089 :
1090 546917 : DO i = 1, nvalues
1091 515376 : value = my_index(i)
1092 515376 : CALL get_indexes(value, ntype, tmpij)
1093 515376 : pot => potparm%pot(tmpij(1), tmpij(2))%pot
1094 :
1095 515376 : IF (do_zbl) THEN
1096 48 : atomic_kind => atomic_kind_set(tmpij(1))
1097 48 : CALL get_atomic_kind(atomic_kind, rcov=rcovi, z=zi)
1098 48 : atomic_kind => atomic_kind_set(tmpij(2))
1099 48 : CALL get_atomic_kind(atomic_kind, rcov=rcovj, z=zj)
1100 : CALL zbl_matching_polinomial(pot, rcovi, rcovj, REAL(zi, KIND=dp), &
1101 48 : REAL(zj, KIND=dp))
1102 : END IF
1103 : ! Derivative factors
1104 1030752 : pot%spl_f%dscale = pot%spl_f%fscale/pot%spl_f%rscale
1105 : ! Cutoff for the potentials on splines
1106 546917 : IF (shift_cutoff) THEN
1107 : ! Cutoff NonBonded
1108 114956 : pot%spl_f%cutoff = ener_pot(pot, hicut0, 0.0_dp)
1109 : END IF
1110 : END DO
1111 :
1112 : ! Handle the cutoff
1113 31541 : IF (shift_cutoff) THEN
1114 28021 : pot_ref => potparm%pot(tmpij_out(1), tmpij_out(2))%pot
1115 142977 : DO i = 1, nvalues
1116 114956 : value = my_index(i)
1117 114956 : CALL get_indexes(value, ntype, tmpij)
1118 114956 : pot => potparm%pot(tmpij(1), tmpij(2))%pot
1119 114956 : IF (value == min_val) CYCLE
1120 : ! Cutoff NonBonded
1121 142977 : pot%spl_f%cutoff = pot_ref%spl_f%cutoff*pot%spl_f%fscale(1) - pot%spl_f%cutoff
1122 : END DO
1123 : END IF
1124 31541 : CALL finalizef()
1125 :
1126 31541 : CALL timestop(handle)
1127 :
1128 31541 : END SUBROUTINE set_potparm_index
1129 :
1130 : ! **************************************************************************************************
1131 : !> \brief Gives back the indices of the matrix w.r.t. the collective array index
1132 : !> \param Inind ...
1133 : !> \param ndim ...
1134 : !> \param ij ...
1135 : !> \author Teodoro Laino [tlaino] 2006.05
1136 : ! **************************************************************************************************
1137 2668715 : SUBROUTINE get_indexes(Inind, ndim, ij)
1138 : INTEGER, INTENT(IN) :: Inind, ndim
1139 : INTEGER, DIMENSION(2), INTENT(OUT) :: ij
1140 :
1141 : INTEGER :: i, tmp
1142 :
1143 2668715 : tmp = 0
1144 8006145 : ij = HUGE(0)
1145 355421855 : DO i = 1, ndim
1146 355421855 : tmp = tmp + i
1147 355421855 : IF (tmp >= Inind) THEN
1148 2668715 : ij(1) = i
1149 2668715 : ij(2) = Inind - tmp + i
1150 2668715 : EXIT
1151 : END IF
1152 : END DO
1153 2668715 : END SUBROUTINE get_indexes
1154 :
1155 : END MODULE pair_potential
1156 :
|
