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