(git:618e178)
Loading...
Searching...
No Matches
eip_silicon.F
Go to the documentation of this file.
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!> \brief Empirical interatomic potentials for Silicon
10!> \note
11!> Stefan Goedecker's OpenMP implementation of Bazant's EDIP & Lenosky's
12!> empirical interatomic potentials for Silicon.
13!> \par History
14!> 03.2006 initial create [tdk]
15!> \author Thomas D. Kuehne (tkuehne@cp2k.org)
16! **************************************************************************************************
17MODULE eip_silicon
18 USE atomic_kind_list_types, ONLY: atomic_kind_list_type
19 USE atomic_kind_types, ONLY: atomic_kind_type,&
20 get_atomic_kind
21 USE cell_types, ONLY: cell_type,&
22 get_cell
23 USE cp_log_handling, ONLY: cp_get_default_logger,&
24 cp_logger_type
25 USE cp_output_handling, ONLY: cp_p_file,&
26 cp_print_key_finished_output,&
27 cp_print_key_should_output,&
28 cp_print_key_unit_nr
29 USE cp_subsys_types, ONLY: cp_subsys_get,&
30 cp_subsys_type
31 USE distribution_1d_types, ONLY: distribution_1d_type
32 USE eip_environment_types, ONLY: eip_env_get,&
33 eip_environment_type
34 USE input_section_types, ONLY: section_vals_get_subs_vals,&
35 section_vals_type
36 USE kinds, ONLY: dp
37 USE message_passing, ONLY: mp_para_env_type
38 USE particle_types, ONLY: particle_type
39 USE physcon, ONLY: angstrom,&
40 evolt
41
42!$ USE OMP_LIB, ONLY: omp_get_max_threads, omp_get_thread_num, omp_get_num_threads
43#include "./base/base_uses.f90"
44
45 IMPLICIT NONE
46 PRIVATE
47
48 CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'eip_silicon'
49
50 ! *** Public subroutines ***
51 PUBLIC :: eip_bazant, eip_lenosky, eip_stillinger_weber, eip_tersoff
52
53!***
54
55CONTAINS
56
57! **************************************************************************************************
58!> \brief Interface routine of Goedecker's Bazant EDIP to CP2K
59!> \param eip_env ...
60!> \par Literature
61!> http://www-math.mit.edu/~bazant/EDIP
62!> M.Z. Bazant & E. Kaxiras: Modeling of Covalent Bonding in Solids by
63!> Inversion of Cohesive Energy Curves;
64!> Phys. Rev. Lett. 77, 4370 (1996)
65!> M.Z. Bazant, E. Kaxiras and J.F. Justo: Environment-dependent interatomic
66!> potential for bulk silicon;
67!> Phys. Rev. B 56, 8542-8552 (1997)
68!> S. Goedecker: Optimization and parallelization of a force field for silicon
69!> using OpenMP; CPC 148, 1 (2002)
70!> \par History
71!> 03.2006 initial create [tdk]
72!> \author Thomas D. Kuehne (tkuehne@cp2k.org)
73! **************************************************************************************************
74 SUBROUTINE eip_bazant(eip_env)
75 TYPE(eip_environment_type), POINTER :: eip_env
76
77 CHARACTER(len=*), PARAMETER :: routinen = 'eip_bazant'
78
79 INTEGER :: handle, i, iparticle, iparticle_kind, &
80 iparticle_local, iw, natom, &
81 nparticle_kind, nparticle_local
82 REAL(kind=dp) :: ekin, ener, ener_var, mass
83 REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :) :: rxyz
84 REAL(kind=dp), DIMENSION(3) :: abc
85 TYPE(atomic_kind_list_type), POINTER :: atomic_kinds
86 TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
87 TYPE(atomic_kind_type), POINTER :: atomic_kind
88 TYPE(cell_type), POINTER :: cell
89 TYPE(cp_logger_type), POINTER :: logger
90 TYPE(cp_subsys_type), POINTER :: subsys
91 TYPE(distribution_1d_type), POINTER :: local_particles
92 TYPE(mp_para_env_type), POINTER :: para_env
93 TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
94 TYPE(section_vals_type), POINTER :: eip_section
95
96! ------------------------------------------------------------------------
97
98 CALL timeset(routinen, handle)
99
100 NULLIFY (cell, particle_set, eip_section, logger, atomic_kinds, &
101 atomic_kind, local_particles, subsys, atomic_kind_set, para_env)
102
103 ekin = 0.0_dp
104
105 logger => cp_get_default_logger()
106
107 cpassert(ASSOCIATED(eip_env))
108
109 CALL eip_env_get(eip_env=eip_env, cell=cell, particle_set=particle_set, &
110 subsys=subsys, local_particles=local_particles, &
111 atomic_kind_set=atomic_kind_set)
112 CALL get_cell(cell=cell, abc=abc)
113
114 eip_section => section_vals_get_subs_vals(eip_env%force_env_input, "EIP")
115 natom = SIZE(particle_set)
116 !natom = local_particles%n_el(1)
117
118 ALLOCATE (rxyz(3, natom))
119
120 DO i = 1, natom
121 !iparticle = local_particles%list(1)%array(i)
122 rxyz(:, i) = particle_set(i)%r(:)*angstrom
123 END DO
124
125 CALL eip_bazant_silicon(nat=natom, alat=abc*angstrom, rxyz0=rxyz, &
126 fxyz=eip_env%eip_forces, ener=ener, &
127 coord=eip_env%coord_avg, ener_var=ener_var, &
128 coord_var=eip_env%coord_var, count=eip_env%count)
129
130 !CALL get_part_ke(md_env, tbmd_energy%E_kinetic, int_grp=globalenv%para_env)
131 CALL cp_subsys_get(subsys=subsys, atomic_kinds=atomic_kinds)
132
133 nparticle_kind = atomic_kinds%n_els
134
135 DO iparticle_kind = 1, nparticle_kind
136 atomic_kind => atomic_kind_set(iparticle_kind)
137 CALL get_atomic_kind(atomic_kind=atomic_kind, mass=mass)
138 nparticle_local = local_particles%n_el(iparticle_kind)
139 DO iparticle_local = 1, nparticle_local
140 iparticle = local_particles%list(iparticle_kind)%array(iparticle_local)
141 ekin = ekin + 0.5_dp*mass* &
142 (particle_set(iparticle)%v(1)*particle_set(iparticle)%v(1) &
143 + particle_set(iparticle)%v(2)*particle_set(iparticle)%v(2) &
144 + particle_set(iparticle)%v(3)*particle_set(iparticle)%v(3))
145 END DO
146 END DO
147
148 ! sum all contributions to energy over calculated parts on all processors
149 CALL cp_subsys_get(subsys=subsys, para_env=para_env)
150 CALL para_env%sum(ekin)
151 eip_env%eip_kinetic_energy = ekin
152
153 eip_env%eip_potential_energy = ener/evolt
154 eip_env%eip_energy = eip_env%eip_kinetic_energy + eip_env%eip_potential_energy
155 eip_env%eip_energy_var = ener_var/evolt
156
157 DO i = 1, natom
158 particle_set(i)%f(:) = eip_env%eip_forces(:, i)/evolt*angstrom
159 END DO
160
161 DEALLOCATE (rxyz)
162
163 ! Print
164 IF (btest(cp_print_key_should_output(logger%iter_info, &
165 eip_section, "PRINT%ENERGIES"), cp_p_file)) THEN
166 iw = cp_print_key_unit_nr(logger, eip_section, "PRINT%ENERGIES", &
167 extension=".mmLog")
168
169 CALL eip_print_energies(eip_env=eip_env, output_unit=iw)
170 CALL cp_print_key_finished_output(iw, logger, eip_section, &
171 "PRINT%ENERGIES")
172 END IF
173
174 IF (btest(cp_print_key_should_output(logger%iter_info, &
175 eip_section, "PRINT%ENERGIES_VAR"), cp_p_file)) THEN
176 iw = cp_print_key_unit_nr(logger, eip_section, "PRINT%ENERGIES_VAR", &
177 extension=".mmLog")
178
179 CALL eip_print_energy_var(eip_env=eip_env, output_unit=iw)
180 CALL cp_print_key_finished_output(iw, logger, eip_section, &
181 "PRINT%ENERGIES_VAR")
182 END IF
183
184 IF (btest(cp_print_key_should_output(logger%iter_info, &
185 eip_section, "PRINT%FORCES"), cp_p_file)) THEN
186 iw = cp_print_key_unit_nr(logger, eip_section, "PRINT%FORCES", &
187 extension=".mmLog")
188
189 CALL eip_print_forces(eip_env=eip_env, output_unit=iw)
190 CALL cp_print_key_finished_output(iw, logger, eip_section, &
191 "PRINT%FORCES")
192 END IF
193
194 IF (btest(cp_print_key_should_output(logger%iter_info, &
195 eip_section, "PRINT%COORD_AVG"), cp_p_file)) THEN
196 iw = cp_print_key_unit_nr(logger, eip_section, "PRINT%COORD_AVG", &
197 extension=".mmLog")
198
199 CALL eip_print_coord_avg(eip_env=eip_env, output_unit=iw)
200 CALL cp_print_key_finished_output(iw, logger, eip_section, &
201 "PRINT%COORD_AVG")
202 END IF
203
204 IF (btest(cp_print_key_should_output(logger%iter_info, &
205 eip_section, "PRINT%COORD_VAR"), cp_p_file)) THEN
206 iw = cp_print_key_unit_nr(logger, eip_section, "PRINT%COORD_VAR", &
207 extension=".mmLog")
208
209 CALL eip_print_coord_var(eip_env=eip_env, output_unit=iw)
210 CALL cp_print_key_finished_output(iw, logger, eip_section, &
211 "PRINT%COORD_VAR")
212 END IF
213
214 IF (btest(cp_print_key_should_output(logger%iter_info, &
215 eip_section, "PRINT%COUNT"), cp_p_file)) THEN
216 iw = cp_print_key_unit_nr(logger, eip_section, "PRINT%COUNT", &
217 extension=".mmLog")
218
219 CALL eip_print_count(eip_env=eip_env, output_unit=iw)
220 CALL cp_print_key_finished_output(iw, logger, eip_section, &
221 "PRINT%COUNT")
222 END IF
223
224 CALL timestop(handle)
225
226 END SUBROUTINE eip_bazant
227
228! **************************************************************************************************
229!> \brief Interface routine of Goedecker's Lenosky force field to CP2K
230!> \param eip_env ...
231!> \par Literature
232!> T. Lenosky, et. al.: Highly optimized empirical potential model of silicon;
233!> Modelling Simul. Sci. Eng., 8 (2000)
234!> S. Goedecker: Optimization and parallelization of a force field for silicon
235!> using OpenMP; CPC 148, 1 (2002)
236!> \par History
237!> 03.2006 initial create [tdk]
238!> \author Thomas D. Kuehne (tkuehne@cp2k.org)
239! **************************************************************************************************
240 SUBROUTINE eip_lenosky(eip_env)
241 TYPE(eip_environment_type), POINTER :: eip_env
242
243 CHARACTER(len=*), PARAMETER :: routinen = 'eip_lenosky'
244
245 INTEGER :: handle, i, iparticle, iparticle_kind, &
246 iparticle_local, iw, natom, &
247 nparticle_kind, nparticle_local
248 REAL(kind=dp) :: ekin, ener, ener_var, mass
249 REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :) :: rxyz
250 REAL(kind=dp), DIMENSION(3) :: abc
251 TYPE(atomic_kind_list_type), POINTER :: atomic_kinds
252 TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
253 TYPE(atomic_kind_type), POINTER :: atomic_kind
254 TYPE(cell_type), POINTER :: cell
255 TYPE(cp_logger_type), POINTER :: logger
256 TYPE(cp_subsys_type), POINTER :: subsys
257 TYPE(distribution_1d_type), POINTER :: local_particles
258 TYPE(mp_para_env_type), POINTER :: para_env
259 TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
260 TYPE(section_vals_type), POINTER :: eip_section
261
262! ------------------------------------------------------------------------
263
264 CALL timeset(routinen, handle)
265
266 NULLIFY (cell, particle_set, eip_section, logger, atomic_kinds, &
267 atomic_kind, local_particles, subsys, atomic_kind_set, para_env)
268
269 ekin = 0.0_dp
270
271 logger => cp_get_default_logger()
272
273 cpassert(ASSOCIATED(eip_env))
274
275 CALL eip_env_get(eip_env=eip_env, cell=cell, particle_set=particle_set, &
276 subsys=subsys, local_particles=local_particles, &
277 atomic_kind_set=atomic_kind_set)
278 CALL get_cell(cell=cell, abc=abc)
279
280 eip_section => section_vals_get_subs_vals(eip_env%force_env_input, "EIP")
281 natom = SIZE(particle_set)
282 !natom = local_particles%n_el(1)
283
284 ALLOCATE (rxyz(3, natom))
285
286 DO i = 1, natom
287 !iparticle = local_particles%list(1)%array(i)
288 rxyz(:, i) = particle_set(i)%r(:)*angstrom
289 END DO
290
291 CALL eip_lenosky_silicon(nat=natom, alat=abc*angstrom, rxyz0=rxyz, &
292 fxyz=eip_env%eip_forces, ener=ener, &
293 coord=eip_env%coord_avg, ener_var=ener_var, &
294 coord_var=eip_env%coord_var, count=eip_env%count)
295
296 !CALL get_part_ke(md_env, tbmd_energy%E_kinetic, int_grp=globalenv%para_env)
297 CALL cp_subsys_get(subsys=subsys, atomic_kinds=atomic_kinds)
298
299 nparticle_kind = atomic_kinds%n_els
300
301 DO iparticle_kind = 1, nparticle_kind
302 atomic_kind => atomic_kind_set(iparticle_kind)
303 CALL get_atomic_kind(atomic_kind=atomic_kind, mass=mass)
304 nparticle_local = local_particles%n_el(iparticle_kind)
305 DO iparticle_local = 1, nparticle_local
306 iparticle = local_particles%list(iparticle_kind)%array(iparticle_local)
307 ekin = ekin + 0.5_dp*mass* &
308 (particle_set(iparticle)%v(1)*particle_set(iparticle)%v(1) &
309 + particle_set(iparticle)%v(2)*particle_set(iparticle)%v(2) &
310 + particle_set(iparticle)%v(3)*particle_set(iparticle)%v(3))
311 END DO
312 END DO
313
314 ! sum all contributions to energy over calculated parts on all processors
315 CALL cp_subsys_get(subsys=subsys, para_env=para_env)
316 CALL para_env%sum(ekin)
317 eip_env%eip_kinetic_energy = ekin
318
319 eip_env%eip_potential_energy = ener/evolt
320 eip_env%eip_energy = eip_env%eip_kinetic_energy + eip_env%eip_potential_energy
321 eip_env%eip_energy_var = ener_var/evolt
322
323 DO i = 1, natom
324 particle_set(i)%f(:) = eip_env%eip_forces(:, i)/evolt*angstrom
325 END DO
326
327 DEALLOCATE (rxyz)
328
329 ! Print
330 IF (btest(cp_print_key_should_output(logger%iter_info, &
331 eip_section, "PRINT%ENERGIES"), cp_p_file)) THEN
332 iw = cp_print_key_unit_nr(logger, eip_section, "PRINT%ENERGIES", &
333 extension=".mmLog")
334
335 CALL eip_print_energies(eip_env=eip_env, output_unit=iw)
336 CALL cp_print_key_finished_output(iw, logger, eip_section, &
337 "PRINT%ENERGIES")
338 END IF
339
340 IF (btest(cp_print_key_should_output(logger%iter_info, &
341 eip_section, "PRINT%ENERGIES_VAR"), cp_p_file)) THEN
342 iw = cp_print_key_unit_nr(logger, eip_section, "PRINT%ENERGIES_VAR", &
343 extension=".mmLog")
344
345 CALL eip_print_energy_var(eip_env=eip_env, output_unit=iw)
346 CALL cp_print_key_finished_output(iw, logger, eip_section, &
347 "PRINT%ENERGIES_VAR")
348 END IF
349
350 IF (btest(cp_print_key_should_output(logger%iter_info, &
351 eip_section, "PRINT%FORCES"), cp_p_file)) THEN
352 iw = cp_print_key_unit_nr(logger, eip_section, "PRINT%FORCES", &
353 extension=".mmLog")
354
355 CALL eip_print_forces(eip_env=eip_env, output_unit=iw)
356 CALL cp_print_key_finished_output(iw, logger, eip_section, &
357 "PRINT%FORCES")
358 END IF
359
360 IF (btest(cp_print_key_should_output(logger%iter_info, &
361 eip_section, "PRINT%COORD_AVG"), cp_p_file)) THEN
362 iw = cp_print_key_unit_nr(logger, eip_section, "PRINT%COORD_AVG", &
363 extension=".mmLog")
364
365 CALL eip_print_coord_avg(eip_env=eip_env, output_unit=iw)
366 CALL cp_print_key_finished_output(iw, logger, eip_section, &
367 "PRINT%COORD_AVG")
368 END IF
369
370 IF (btest(cp_print_key_should_output(logger%iter_info, &
371 eip_section, "PRINT%COORD_VAR"), cp_p_file)) THEN
372 iw = cp_print_key_unit_nr(logger, eip_section, "PRINT%COORD_VAR", &
373 extension=".mmLog")
374
375 CALL eip_print_coord_var(eip_env=eip_env, output_unit=iw)
376 CALL cp_print_key_finished_output(iw, logger, eip_section, &
377 "PRINT%COORD_VAR")
378 END IF
379
380 IF (btest(cp_print_key_should_output(logger%iter_info, &
381 eip_section, "PRINT%COUNT"), cp_p_file)) THEN
382 iw = cp_print_key_unit_nr(logger, eip_section, "PRINT%COUNT", &
383 extension=".mmLog")
384
385 CALL eip_print_count(eip_env=eip_env, output_unit=iw)
386 CALL cp_print_key_finished_output(iw, logger, eip_section, &
387 "PRINT%COUNT")
388 END IF
389
390 CALL timestop(handle)
391
392 END SUBROUTINE eip_lenosky
393
394! **************************************************************************************************
395!> \brief Interface routine of the Stillinger-Weber force field to CP2K
396!> \param eip_env ...
397!> \par Literature
398!> F.H. Stillinger and T.A. Weber:
399!> Computer simulation of local order in condensed phases of silicon;
400!> Phys. Rev. B 31, 5262 (1985)
401!> \par History
402!> 04.2026 added [Thomas D. Kuehne, tkuehne@cp2k.org]
403! **************************************************************************************************
404 SUBROUTINE eip_stillinger_weber(eip_env)
405 TYPE(eip_environment_type), POINTER :: eip_env
406
407 CHARACTER(len=*), PARAMETER :: routinen = 'eip_stillinger_weber'
408
409 INTEGER :: handle, i, iparticle, iparticle_kind, &
410 iparticle_local, iw, natom, &
411 nparticle_kind, nparticle_local
412 REAL(kind=dp) :: ekin, ener, mass
413 REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :) :: rxyz
414 REAL(kind=dp), DIMENSION(3) :: abc
415 TYPE(atomic_kind_list_type), POINTER :: atomic_kinds
416 TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
417 TYPE(atomic_kind_type), POINTER :: atomic_kind
418 TYPE(cell_type), POINTER :: cell
419 TYPE(cp_logger_type), POINTER :: logger
420 TYPE(cp_subsys_type), POINTER :: subsys
421 TYPE(distribution_1d_type), POINTER :: local_particles
422 TYPE(mp_para_env_type), POINTER :: para_env
423 TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
424 TYPE(section_vals_type), POINTER :: eip_section
425
426 CALL timeset(routinen, handle)
427
428 NULLIFY (cell, particle_set, eip_section, logger, atomic_kinds, &
429 atomic_kind, local_particles, subsys, atomic_kind_set, para_env)
430
431 ekin = 0.0_dp
432
433 logger => cp_get_default_logger()
434
435 cpassert(ASSOCIATED(eip_env))
436
437 CALL eip_env_get(eip_env=eip_env, cell=cell, particle_set=particle_set, &
438 subsys=subsys, local_particles=local_particles, &
439 atomic_kind_set=atomic_kind_set)
440 CALL get_cell(cell=cell, abc=abc)
441
442 eip_section => section_vals_get_subs_vals(eip_env%force_env_input, "EIP")
443 natom = SIZE(particle_set)
444
445 ALLOCATE (rxyz(3, natom))
446
447 DO i = 1, natom
448 rxyz(:, i) = particle_set(i)%r(:)*angstrom
449 END DO
450
451 CALL eip_stillinger_weber_silicon(nat=natom, alat=abc*angstrom, &
452 rxyz0=rxyz, fxyz=eip_env%eip_forces, &
453 etot=ener, count=eip_env%count)
454
455 eip_env%coord_avg = 0.0_dp
456 eip_env%coord_var = 0.0_dp
457
458 CALL cp_subsys_get(subsys=subsys, atomic_kinds=atomic_kinds)
459
460 nparticle_kind = atomic_kinds%n_els
461
462 DO iparticle_kind = 1, nparticle_kind
463 atomic_kind => atomic_kind_set(iparticle_kind)
464 CALL get_atomic_kind(atomic_kind=atomic_kind, mass=mass)
465 nparticle_local = local_particles%n_el(iparticle_kind)
466 DO iparticle_local = 1, nparticle_local
467 iparticle = local_particles%list(iparticle_kind)%array(iparticle_local)
468 ekin = ekin + 0.5_dp*mass* &
469 (particle_set(iparticle)%v(1)*particle_set(iparticle)%v(1) &
470 + particle_set(iparticle)%v(2)*particle_set(iparticle)%v(2) &
471 + particle_set(iparticle)%v(3)*particle_set(iparticle)%v(3))
472 END DO
473 END DO
474
475 CALL cp_subsys_get(subsys=subsys, para_env=para_env)
476 CALL para_env%sum(ekin)
477 eip_env%eip_kinetic_energy = ekin
478
479 eip_env%eip_potential_energy = ener/evolt
480 eip_env%eip_energy = eip_env%eip_kinetic_energy + eip_env%eip_potential_energy
481 eip_env%eip_energy_var = 0.0_dp
482
483 DO i = 1, natom
484 particle_set(i)%f(:) = eip_env%eip_forces(:, i)/evolt*angstrom
485 END DO
486
487 DEALLOCATE (rxyz)
488
489 IF (btest(cp_print_key_should_output(logger%iter_info, &
490 eip_section, "PRINT%ENERGIES"), cp_p_file)) THEN
491 iw = cp_print_key_unit_nr(logger, eip_section, "PRINT%ENERGIES", &
492 extension=".mmLog")
493
494 CALL eip_print_energies(eip_env=eip_env, output_unit=iw)
495 CALL cp_print_key_finished_output(iw, logger, eip_section, &
496 "PRINT%ENERGIES")
497 END IF
498
499 IF (btest(cp_print_key_should_output(logger%iter_info, &
500 eip_section, "PRINT%ENERGIES_VAR"), cp_p_file)) THEN
501 iw = cp_print_key_unit_nr(logger, eip_section, "PRINT%ENERGIES_VAR", &
502 extension=".mmLog")
503
504 CALL eip_print_energy_var(eip_env=eip_env, output_unit=iw)
505 CALL cp_print_key_finished_output(iw, logger, eip_section, &
506 "PRINT%ENERGIES_VAR")
507 END IF
508
509 IF (btest(cp_print_key_should_output(logger%iter_info, &
510 eip_section, "PRINT%FORCES"), cp_p_file)) THEN
511 iw = cp_print_key_unit_nr(logger, eip_section, "PRINT%FORCES", &
512 extension=".mmLog")
513
514 CALL eip_print_forces(eip_env=eip_env, output_unit=iw)
515 CALL cp_print_key_finished_output(iw, logger, eip_section, &
516 "PRINT%FORCES")
517 END IF
518
519 IF (btest(cp_print_key_should_output(logger%iter_info, &
520 eip_section, "PRINT%COORD_AVG"), cp_p_file)) THEN
521 iw = cp_print_key_unit_nr(logger, eip_section, "PRINT%COORD_AVG", &
522 extension=".mmLog")
523
524 CALL eip_print_coord_avg(eip_env=eip_env, output_unit=iw)
525 CALL cp_print_key_finished_output(iw, logger, eip_section, &
526 "PRINT%COORD_AVG")
527 END IF
528
529 IF (btest(cp_print_key_should_output(logger%iter_info, &
530 eip_section, "PRINT%COORD_VAR"), cp_p_file)) THEN
531 iw = cp_print_key_unit_nr(logger, eip_section, "PRINT%COORD_VAR", &
532 extension=".mmLog")
533
534 CALL eip_print_coord_var(eip_env=eip_env, output_unit=iw)
535 CALL cp_print_key_finished_output(iw, logger, eip_section, &
536 "PRINT%COORD_VAR")
537 END IF
538
539 IF (btest(cp_print_key_should_output(logger%iter_info, &
540 eip_section, "PRINT%COUNT"), cp_p_file)) THEN
541 iw = cp_print_key_unit_nr(logger, eip_section, "PRINT%COUNT", &
542 extension=".mmLog")
543
544 CALL eip_print_count(eip_env=eip_env, output_unit=iw)
545 CALL cp_print_key_finished_output(iw, logger, eip_section, &
546 "PRINT%COUNT")
547 END IF
548
549 CALL timestop(handle)
550
551 END SUBROUTINE eip_stillinger_weber
552
553! **************************************************************************************************
554!> \brief Interface routine of the Tersoff force field to CP2K
555!> \param eip_env ...
556!> \par Literature
557!> J. Tersoff:
558!> New empirical approach for the structure and energy of covalent systems;
559!> Phys. Rev. Lett. 61, 2879 (1988)
560!> J. Tersoff:
561!> Modeling solid-state chemistry: Interatomic potentials for multicomponent systems;
562!> Phys. Rev. B 39, 5566 (1989)
563!> \par History
564!> 04.2026 added [Thomas D. Kuehne, tkuehne@cp2k.org]
565! **************************************************************************************************
566 SUBROUTINE eip_tersoff(eip_env)
567 TYPE(eip_environment_type), POINTER :: eip_env
568
569 CHARACTER(len=*), PARAMETER :: routinen = 'eip_tersoff'
570
571 INTEGER :: handle, i, iparticle, iparticle_kind, &
572 iparticle_local, iw, natom, &
573 nparticle_kind, nparticle_local
574 REAL(kind=dp) :: ekin, ener, mass
575 REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :) :: rxyz
576 REAL(kind=dp), DIMENSION(3) :: abc
577 TYPE(atomic_kind_list_type), POINTER :: atomic_kinds
578 TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
579 TYPE(atomic_kind_type), POINTER :: atomic_kind
580 TYPE(cell_type), POINTER :: cell
581 TYPE(cp_logger_type), POINTER :: logger
582 TYPE(cp_subsys_type), POINTER :: subsys
583 TYPE(distribution_1d_type), POINTER :: local_particles
584 TYPE(mp_para_env_type), POINTER :: para_env
585 TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
586 TYPE(section_vals_type), POINTER :: eip_section
587
588 CALL timeset(routinen, handle)
589
590 NULLIFY (cell, particle_set, eip_section, logger, atomic_kinds, &
591 atomic_kind, local_particles, subsys, atomic_kind_set, para_env)
592
593 ekin = 0.0_dp
594
595 logger => cp_get_default_logger()
596
597 cpassert(ASSOCIATED(eip_env))
598
599 CALL eip_env_get(eip_env=eip_env, cell=cell, particle_set=particle_set, &
600 subsys=subsys, local_particles=local_particles, &
601 atomic_kind_set=atomic_kind_set)
602 CALL get_cell(cell=cell, abc=abc)
603
604 eip_section => section_vals_get_subs_vals(eip_env%force_env_input, "EIP")
605 natom = SIZE(particle_set)
606
607 ALLOCATE (rxyz(3, natom))
608
609 DO i = 1, natom
610 rxyz(:, i) = particle_set(i)%r(:)*angstrom
611 END DO
612
613 CALL eip_tersoff_silicon(nat=natom, alat=abc*angstrom, rxyz=rxyz, &
614 fxyz=eip_env%eip_forces, etot=ener, &
615 count=eip_env%count)
616
617 eip_env%coord_avg = 0.0_dp
618 eip_env%coord_var = 0.0_dp
619
620 CALL cp_subsys_get(subsys=subsys, atomic_kinds=atomic_kinds)
621
622 nparticle_kind = atomic_kinds%n_els
623
624 DO iparticle_kind = 1, nparticle_kind
625 atomic_kind => atomic_kind_set(iparticle_kind)
626 CALL get_atomic_kind(atomic_kind=atomic_kind, mass=mass)
627 nparticle_local = local_particles%n_el(iparticle_kind)
628 DO iparticle_local = 1, nparticle_local
629 iparticle = local_particles%list(iparticle_kind)%array(iparticle_local)
630 ekin = ekin + 0.5_dp*mass* &
631 (particle_set(iparticle)%v(1)*particle_set(iparticle)%v(1) &
632 + particle_set(iparticle)%v(2)*particle_set(iparticle)%v(2) &
633 + particle_set(iparticle)%v(3)*particle_set(iparticle)%v(3))
634 END DO
635 END DO
636
637 CALL cp_subsys_get(subsys=subsys, para_env=para_env)
638 CALL para_env%sum(ekin)
639 eip_env%eip_kinetic_energy = ekin
640
641 eip_env%eip_potential_energy = ener/evolt
642 eip_env%eip_energy = eip_env%eip_kinetic_energy + eip_env%eip_potential_energy
643 eip_env%eip_energy_var = 0.0_dp
644
645 DO i = 1, natom
646 particle_set(i)%f(:) = eip_env%eip_forces(:, i)/evolt*angstrom
647 END DO
648
649 DEALLOCATE (rxyz)
650
651 IF (btest(cp_print_key_should_output(logger%iter_info, &
652 eip_section, "PRINT%ENERGIES"), cp_p_file)) THEN
653 iw = cp_print_key_unit_nr(logger, eip_section, "PRINT%ENERGIES", &
654 extension=".mmLog")
655
656 CALL eip_print_energies(eip_env=eip_env, output_unit=iw)
657 CALL cp_print_key_finished_output(iw, logger, eip_section, &
658 "PRINT%ENERGIES")
659 END IF
660
661 IF (btest(cp_print_key_should_output(logger%iter_info, &
662 eip_section, "PRINT%ENERGIES_VAR"), cp_p_file)) THEN
663 iw = cp_print_key_unit_nr(logger, eip_section, "PRINT%ENERGIES_VAR", &
664 extension=".mmLog")
665
666 CALL eip_print_energy_var(eip_env=eip_env, output_unit=iw)
667 CALL cp_print_key_finished_output(iw, logger, eip_section, &
668 "PRINT%ENERGIES_VAR")
669 END IF
670
671 IF (btest(cp_print_key_should_output(logger%iter_info, &
672 eip_section, "PRINT%FORCES"), cp_p_file)) THEN
673 iw = cp_print_key_unit_nr(logger, eip_section, "PRINT%FORCES", &
674 extension=".mmLog")
675
676 CALL eip_print_forces(eip_env=eip_env, output_unit=iw)
677 CALL cp_print_key_finished_output(iw, logger, eip_section, &
678 "PRINT%FORCES")
679 END IF
680
681 IF (btest(cp_print_key_should_output(logger%iter_info, &
682 eip_section, "PRINT%COORD_AVG"), cp_p_file)) THEN
683 iw = cp_print_key_unit_nr(logger, eip_section, "PRINT%COORD_AVG", &
684 extension=".mmLog")
685
686 CALL eip_print_coord_avg(eip_env=eip_env, output_unit=iw)
687 CALL cp_print_key_finished_output(iw, logger, eip_section, &
688 "PRINT%COORD_AVG")
689 END IF
690
691 IF (btest(cp_print_key_should_output(logger%iter_info, &
692 eip_section, "PRINT%COORD_VAR"), cp_p_file)) THEN
693 iw = cp_print_key_unit_nr(logger, eip_section, "PRINT%COORD_VAR", &
694 extension=".mmLog")
695
696 CALL eip_print_coord_var(eip_env=eip_env, output_unit=iw)
697 CALL cp_print_key_finished_output(iw, logger, eip_section, &
698 "PRINT%COORD_VAR")
699 END IF
700
701 IF (btest(cp_print_key_should_output(logger%iter_info, &
702 eip_section, "PRINT%COUNT"), cp_p_file)) THEN
703 iw = cp_print_key_unit_nr(logger, eip_section, "PRINT%COUNT", &
704 extension=".mmLog")
705
706 CALL eip_print_count(eip_env=eip_env, output_unit=iw)
707 CALL cp_print_key_finished_output(iw, logger, eip_section, &
708 "PRINT%COUNT")
709 END IF
710
711 CALL timestop(handle)
712
713 END SUBROUTINE eip_tersoff
714
715! **************************************************************************************************
716!> \brief Print routine for the EIP energies
717!> \param eip_env The eip environment of matter
718!> \param output_unit The output unit
719!> \par History
720!> 03.2006 initial create [tdk]
721!> \author Thomas D. Kuehne (tkuehne@cp2k.org)
722!> \note
723!> As usual the EIP energies differ from the DFT energies!
724!> Only the relative energy differences are correctly reproduced.
725! **************************************************************************************************
726 SUBROUTINE eip_print_energies(eip_env, output_unit)
727 TYPE(eip_environment_type), POINTER :: eip_env
728 INTEGER, INTENT(IN) :: output_unit
729
730! ------------------------------------------------------------------------
731
732 IF (output_unit > 0) THEN
733 WRITE (unit=output_unit, fmt="(/,(T3,A,T55,F25.14))") &
734 "Kinetic energy [Hartree]: ", eip_env%eip_kinetic_energy, &
735 "Potential energy [Hartree]: ", eip_env%eip_potential_energy, &
736 "Total EIP energy [Hartree]: ", eip_env%eip_energy
737 END IF
738
739 END SUBROUTINE eip_print_energies
740
741! **************************************************************************************************
742!> \brief Print routine for the variance of the energy/atom
743!> \param eip_env The eip environment of matter
744!> \param output_unit The output unit
745!> \par History
746!> 03.2006 initial create [tdk]
747!> \author Thomas D. Kuehne (tkuehne@cp2k.org)
748! **************************************************************************************************
749 SUBROUTINE eip_print_energy_var(eip_env, output_unit)
750 TYPE(eip_environment_type), POINTER :: eip_env
751 INTEGER, INTENT(IN) :: output_unit
752
753 INTEGER :: unit_nr
754
755! ------------------------------------------------------------------------
756
757 unit_nr = output_unit
758
759 IF (unit_nr > 0) THEN
760
761 WRITE (unit_nr, *) ""
762 WRITE (unit_nr, *) "The variance of the EIP energy/atom!"
763 WRITE (unit_nr, *) ""
764 WRITE (unit_nr, *) eip_env%eip_energy_var
765
766 END IF
767
768 END SUBROUTINE eip_print_energy_var
769
770! **************************************************************************************************
771!> \brief Print routine for the forces
772!> \param eip_env The eip environment of matter
773!> \param output_unit The output unit
774!> \par History
775!> 03.2006 initial create [tdk]
776!> \author Thomas D. Kuehne (tkuehne@cp2k.org)
777! **************************************************************************************************
778 SUBROUTINE eip_print_forces(eip_env, output_unit)
779 TYPE(eip_environment_type), POINTER :: eip_env
780 INTEGER, INTENT(IN) :: output_unit
781
782 INTEGER :: iatom, natom, unit_nr
783 TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
784
785! ------------------------------------------------------------------------
786
787 NULLIFY (particle_set)
788
789 unit_nr = output_unit
790
791 IF (unit_nr > 0) THEN
792
793 CALL eip_env_get(eip_env=eip_env, particle_set=particle_set)
794
795 natom = SIZE(particle_set)
796
797 WRITE (unit_nr, *) ""
798 WRITE (unit_nr, *) "The EIP forces!"
799 WRITE (unit_nr, *) ""
800 WRITE (unit_nr, *) "Total EIP forces [Hartree/Bohr]"
801 DO iatom = 1, natom
802 WRITE (unit_nr, *) eip_env%eip_forces(1:3, iatom)
803 END DO
804
805 END IF
806
807 END SUBROUTINE eip_print_forces
808
809! **************************************************************************************************
810!> \brief Print routine for the average coordination number
811!> \param eip_env The eip environment of matter
812!> \param output_unit The output unit
813!> \par History
814!> 03.2006 initial create [tdk]
815!> \author Thomas D. Kuehne (tkuehne@cp2k.org)
816! **************************************************************************************************
817 SUBROUTINE eip_print_coord_avg(eip_env, output_unit)
818 TYPE(eip_environment_type), POINTER :: eip_env
819 INTEGER, INTENT(IN) :: output_unit
820
821 INTEGER :: unit_nr
822
823! ------------------------------------------------------------------------
824
825 unit_nr = output_unit
826
827 IF (unit_nr > 0) THEN
828
829 WRITE (unit_nr, *) ""
830 WRITE (unit_nr, *) "The average coordination number!"
831 WRITE (unit_nr, *) ""
832 WRITE (unit_nr, *) eip_env%coord_avg
833
834 END IF
835
836 END SUBROUTINE eip_print_coord_avg
837
838! **************************************************************************************************
839!> \brief Print routine for the variance of the coordination number
840!> \param eip_env The eip environment of matter
841!> \param output_unit The output unit
842!> \par History
843!> 03.2006 initial create [tdk]
844!> \author Thomas D. Kuehne (tkuehne@cp2k.org)
845! **************************************************************************************************
846 SUBROUTINE eip_print_coord_var(eip_env, output_unit)
847 TYPE(eip_environment_type), POINTER :: eip_env
848 INTEGER, INTENT(IN) :: output_unit
849
850 INTEGER :: unit_nr
851
852! ------------------------------------------------------------------------
853
854 unit_nr = output_unit
855
856 IF (unit_nr > 0) THEN
857
858 WRITE (unit_nr, *) ""
859 WRITE (unit_nr, *) "The variance of the coordination number!"
860 WRITE (unit_nr, *) ""
861 WRITE (unit_nr, *) eip_env%coord_var
862
863 END IF
864
865 END SUBROUTINE eip_print_coord_var
866
867! **************************************************************************************************
868!> \brief Print routine for the function call counter
869!> \param eip_env The eip environment of matter
870!> \param output_unit The output unit
871!> \par History
872!> 03.2006 initial create [tdk]
873!> \author Thomas D. Kuehne (tkuehne@cp2k.org)
874! **************************************************************************************************
875 SUBROUTINE eip_print_count(eip_env, output_unit)
876 TYPE(eip_environment_type), POINTER :: eip_env
877 INTEGER, INTENT(IN) :: output_unit
878
879 INTEGER :: unit_nr
880
881! ------------------------------------------------------------------------
882
883 unit_nr = output_unit
884
885 IF (unit_nr > 0) THEN
886
887 WRITE (unit_nr, *) ""
888 WRITE (unit_nr, *) "The function call counter!"
889 WRITE (unit_nr, *) ""
890 WRITE (unit_nr, *) eip_env%count
891
892 END IF
893
894 END SUBROUTINE eip_print_count
895
896! **************************************************************************************************
897!> \brief Bazant's EDIP (environment-dependent interatomic potential) for Silicon
898!> by Stefan Goedecker
899!> \param nat number of atoms
900!> \param alat lattice constants of the orthorombic box containing the particles
901!> \param rxyz0 atomic positions in Angstrom, may be modified on output.
902!> If an atom is outside the box the program will bring it back
903!> into the box by translations through alat
904!> \param fxyz forces in eV/A
905!> \param ener total energy in eV
906!> \param coord average coordination number
907!> \param ener_var variance of the energy/atom
908!> \param coord_var variance of the coordination number
909!> \param count count is increased by one per call, has to be initialized
910!> to 0.e0_dp before first call of eip_bazant
911!> \par Literature
912!> http://www-math.mit.edu/~bazant/EDIP
913!> M.Z. Bazant & E. Kaxiras: Modeling of Covalent Bonding in Solids by
914!> Inversion of Cohesive Energy Curves;
915!> Phys. Rev. Lett. 77, 4370 (1996)
916!> M.Z. Bazant, E. Kaxiras and J.F. Justo: Environment-dependent interatomic
917!> potential for bulk silicon;
918!> Phys. Rev. B 56, 8542-8552 (1997)
919!> S. Goedecker: Optimization and parallelization of a force field for silicon
920!> using OpenMP; CPC 148, 1 (2002)
921!> \par History
922!> 03.2006 initial create [tdk]
923!> \author Thomas D. Kuehne (tkuehne@cp2k.org)
924! **************************************************************************************************
925 SUBROUTINE eip_bazant_silicon(nat, alat, rxyz0, fxyz, ener, coord, ener_var, &
926 coord_var, count)
927
928 INTEGER :: nat
929 REAL(kind=dp) :: alat, rxyz0, fxyz, ener, coord, &
930 ener_var, coord_var, count
931
932 dimension rxyz0(3, nat), fxyz(3, nat), alat(3)
933 REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :) :: rxyz
934 INTEGER, ALLOCATABLE, DIMENSION(:, :) :: lsta
935 INTEGER, ALLOCATABLE, DIMENSION(:) :: lstb
936 INTEGER, ALLOCATABLE, DIMENSION(:) :: lay
937 INTEGER, ALLOCATABLE, DIMENSION(:, :, :, :) :: icell
938 REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :) :: rel
939 REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :) :: txyz
940 REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :) :: s2, s3, sz
941 INTEGER, ALLOCATABLE, DIMENSION(:) :: num2, num3, numz
942
943 REAL(kind=dp) :: coord2, cut, cut2, ener2, rlc1i, rlc2i, rlc3i, tcoord, &
944 tcoord2, tener, tener2
945 INTEGER :: iam, iat, iat1, iat2, ii, i, il, in, indlst, indlstx, istop, &
946 istopg, l2, l3, laymx, ll1, ll2, ll3, lot, max_nbrs, myspace, &
947 l1, myspaceout, ncx, nn, nnbrx, npr
948
949! cut=par_a
950 cut = 3.1213820e0_dp + 1.e-14_dp
951
952 IF (count == 0) OPEN (unit=10, file='bazant.mon', status='unknown')
953 count = count + 1.e0_dp
954
955! linear scaling calculation of verlet list
956 ll1 = int(alat(1)/cut)
957 IF (ll1 < 1) cpabort("alat(1) too small")
958 ll2 = int(alat(2)/cut)
959 IF (ll2 < 1) cpabort("alat(2) too small")
960 ll3 = int(alat(3)/cut)
961 IF (ll3 < 1) cpabort("alat(3) too small")
962
963! determine number of threads
964 npr = 1
965!$OMP PARALLEL PRIVATE(iam) SHARED (npr) DEFAULT(NONE)
966!$ iam = omp_get_thread_num()
967!$ if (iam .eq. 0) npr = omp_get_num_threads()
968!$OMP END PARALLEL
969
970! linear scaling calculation of verlet list
971
972 IF (npr <= 1) THEN !serial if too few processors to gain by parallelizing
973
974! set ncx for serial case, ncx for parallel case set below
975 ncx = 16
976 loop_ncx_s: DO
977 ALLOCATE (icell(0:ncx, -1:ll1, -1:ll2, -1:ll3))
978 icell(0, -1:ll1, -1:ll2, -1:ll3) = 0
979 rlc1i = ll1/alat(1)
980 rlc2i = ll2/alat(2)
981 rlc3i = ll3/alat(3)
982
983 loop_iat_s: DO iat = 1, nat
984 rxyz0(1, iat) = modulo(modulo(rxyz0(1, iat), alat(1)), alat(1))
985 rxyz0(2, iat) = modulo(modulo(rxyz0(2, iat), alat(2)), alat(2))
986 rxyz0(3, iat) = modulo(modulo(rxyz0(3, iat), alat(3)), alat(3))
987 l1 = int(rxyz0(1, iat)*rlc1i)
988 l2 = int(rxyz0(2, iat)*rlc2i)
989 l3 = int(rxyz0(3, iat)*rlc3i)
990
991 ii = icell(0, l1, l2, l3)
992 ii = ii + 1
993 icell(0, l1, l2, l3) = ii
994 IF (ii > ncx) THEN
995 WRITE (10, *) count, 'NCX too small', ncx
996 DEALLOCATE (icell)
997 ncx = ncx*2
998 cycle loop_ncx_s
999 END IF
1000 icell(ii, l1, l2, l3) = iat
1001 END DO loop_iat_s
1002 EXIT loop_ncx_s
1003 END DO loop_ncx_s
1004
1005 ELSE ! parallel case
1006
1007! periodization of particles can be done in parallel
1008!$OMP PARALLEL DO SHARED (alat,nat,rxyz0) PRIVATE(iat) DEFAULT(NONE)
1009 DO iat = 1, nat
1010 rxyz0(1, iat) = modulo(modulo(rxyz0(1, iat), alat(1)), alat(1))
1011 rxyz0(2, iat) = modulo(modulo(rxyz0(2, iat), alat(2)), alat(2))
1012 rxyz0(3, iat) = modulo(modulo(rxyz0(3, iat), alat(3)), alat(3))
1013 END DO
1014!$OMP END PARALLEL DO
1015
1016! assignment to cell is done serially
1017! set ncx for parallel case, ncx for serial case set above
1018 ncx = 16
1019 loop_ncx_p: DO
1020 ALLOCATE (icell(0:ncx, -1:ll1, -1:ll2, -1:ll3))
1021 icell(0, -1:ll1, -1:ll2, -1:ll3) = 0
1022
1023 rlc1i = ll1/alat(1)
1024 rlc2i = ll2/alat(2)
1025 rlc3i = ll3/alat(3)
1026
1027 loop_iat_p: DO iat = 1, nat
1028 l1 = int(rxyz0(1, iat)*rlc1i)
1029 l2 = int(rxyz0(2, iat)*rlc2i)
1030 l3 = int(rxyz0(3, iat)*rlc3i)
1031 ii = icell(0, l1, l2, l3)
1032 ii = ii + 1
1033 icell(0, l1, l2, l3) = ii
1034 IF (ii > ncx) THEN
1035 WRITE (10, *) count, 'NCX too small', ncx
1036 DEALLOCATE (icell)
1037 ncx = ncx*2
1038 cycle loop_ncx_p
1039 END IF
1040 icell(ii, l1, l2, l3) = iat
1041 END DO loop_iat_p
1042 EXIT loop_ncx_p
1043 END DO loop_ncx_p
1044
1045 END IF
1046
1047! duplicate all atoms within boundary layer
1048 laymx = ncx*(2*ll1*ll2 + 2*ll1*ll3 + 2*ll2*ll3 + 4*ll1 + 4*ll2 + 4*ll3 + 8)
1049 nn = nat + laymx
1050 ALLOCATE (rxyz(3, nn), lay(nn))
1051 DO iat = 1, nat
1052 lay(iat) = iat
1053 rxyz(1, iat) = rxyz0(1, iat)
1054 rxyz(2, iat) = rxyz0(2, iat)
1055 rxyz(3, iat) = rxyz0(3, iat)
1056 END DO
1057 il = nat
1058! xy plane
1059 DO l2 = 0, ll2 - 1
1060 DO l1 = 0, ll1 - 1
1061
1062 in = icell(0, l1, l2, 0)
1063 icell(0, l1, l2, ll3) = in
1064 DO ii = 1, in
1065 i = icell(ii, l1, l2, 0)
1066 il = il + 1
1067 IF (il > nn) cpabort("enlarge laymx")
1068 lay(il) = i
1069 icell(ii, l1, l2, ll3) = il
1070 rxyz(1, il) = rxyz(1, i)
1071 rxyz(2, il) = rxyz(2, i)
1072 rxyz(3, il) = rxyz(3, i) + alat(3)
1073 END DO
1074
1075 in = icell(0, l1, l2, ll3 - 1)
1076 icell(0, l1, l2, -1) = in
1077 DO ii = 1, in
1078 i = icell(ii, l1, l2, ll3 - 1)
1079 il = il + 1
1080 IF (il > nn) cpabort("enlarge laymx")
1081 lay(il) = i
1082 icell(ii, l1, l2, -1) = il
1083 rxyz(1, il) = rxyz(1, i)
1084 rxyz(2, il) = rxyz(2, i)
1085 rxyz(3, il) = rxyz(3, i) - alat(3)
1086 END DO
1087
1088 END DO
1089 END DO
1090
1091! yz plane
1092 DO l3 = 0, ll3 - 1
1093 DO l2 = 0, ll2 - 1
1094
1095 in = icell(0, 0, l2, l3)
1096 icell(0, ll1, l2, l3) = in
1097 DO ii = 1, in
1098 i = icell(ii, 0, l2, l3)
1099 il = il + 1
1100 IF (il > nn) cpabort("enlarge laymx")
1101 lay(il) = i
1102 icell(ii, ll1, l2, l3) = il
1103 rxyz(1, il) = rxyz(1, i) + alat(1)
1104 rxyz(2, il) = rxyz(2, i)
1105 rxyz(3, il) = rxyz(3, i)
1106 END DO
1107
1108 in = icell(0, ll1 - 1, l2, l3)
1109 icell(0, -1, l2, l3) = in
1110 DO ii = 1, in
1111 i = icell(ii, ll1 - 1, l2, l3)
1112 il = il + 1
1113 IF (il > nn) cpabort("enlarge laymx")
1114 lay(il) = i
1115 icell(ii, -1, l2, l3) = il
1116 rxyz(1, il) = rxyz(1, i) - alat(1)
1117 rxyz(2, il) = rxyz(2, i)
1118 rxyz(3, il) = rxyz(3, i)
1119 END DO
1120
1121 END DO
1122 END DO
1123
1124! xz plane
1125 DO l3 = 0, ll3 - 1
1126 DO l1 = 0, ll1 - 1
1127
1128 in = icell(0, l1, 0, l3)
1129 icell(0, l1, ll2, l3) = in
1130 DO ii = 1, in
1131 i = icell(ii, l1, 0, l3)
1132 il = il + 1
1133 IF (il > nn) cpabort("enlarge laymx")
1134 lay(il) = i
1135 icell(ii, l1, ll2, l3) = il
1136 rxyz(1, il) = rxyz(1, i)
1137 rxyz(2, il) = rxyz(2, i) + alat(2)
1138 rxyz(3, il) = rxyz(3, i)
1139 END DO
1140
1141 in = icell(0, l1, ll2 - 1, l3)
1142 icell(0, l1, -1, l3) = in
1143 DO ii = 1, in
1144 i = icell(ii, l1, ll2 - 1, l3)
1145 il = il + 1
1146 IF (il > nn) cpabort("enlarge laymx")
1147 lay(il) = i
1148 icell(ii, l1, -1, l3) = il
1149 rxyz(1, il) = rxyz(1, i)
1150 rxyz(2, il) = rxyz(2, i) - alat(2)
1151 rxyz(3, il) = rxyz(3, i)
1152 END DO
1153
1154 END DO
1155 END DO
1156
1157! x axis
1158 DO l1 = 0, ll1 - 1
1159
1160 in = icell(0, l1, 0, 0)
1161 icell(0, l1, ll2, ll3) = in
1162 DO ii = 1, in
1163 i = icell(ii, l1, 0, 0)
1164 il = il + 1
1165 IF (il > nn) cpabort("enlarge laymx")
1166 lay(il) = i
1167 icell(ii, l1, ll2, ll3) = il
1168 rxyz(1, il) = rxyz(1, i)
1169 rxyz(2, il) = rxyz(2, i) + alat(2)
1170 rxyz(3, il) = rxyz(3, i) + alat(3)
1171 END DO
1172
1173 in = icell(0, l1, 0, ll3 - 1)
1174 icell(0, l1, ll2, -1) = in
1175 DO ii = 1, in
1176 i = icell(ii, l1, 0, ll3 - 1)
1177 il = il + 1
1178 IF (il > nn) cpabort("enlarge laymx")
1179 lay(il) = i
1180 icell(ii, l1, ll2, -1) = il
1181 rxyz(1, il) = rxyz(1, i)
1182 rxyz(2, il) = rxyz(2, i) + alat(2)
1183 rxyz(3, il) = rxyz(3, i) - alat(3)
1184 END DO
1185
1186 in = icell(0, l1, ll2 - 1, 0)
1187 icell(0, l1, -1, ll3) = in
1188 DO ii = 1, in
1189 i = icell(ii, l1, ll2 - 1, 0)
1190 il = il + 1
1191 IF (il > nn) cpabort("enlarge laymx")
1192 lay(il) = i
1193 icell(ii, l1, -1, ll3) = il
1194 rxyz(1, il) = rxyz(1, i)
1195 rxyz(2, il) = rxyz(2, i) - alat(2)
1196 rxyz(3, il) = rxyz(3, i) + alat(3)
1197 END DO
1198
1199 in = icell(0, l1, ll2 - 1, ll3 - 1)
1200 icell(0, l1, -1, -1) = in
1201 DO ii = 1, in
1202 i = icell(ii, l1, ll2 - 1, ll3 - 1)
1203 il = il + 1
1204 IF (il > nn) cpabort("enlarge laymx")
1205 lay(il) = i
1206 icell(ii, l1, -1, -1) = il
1207 rxyz(1, il) = rxyz(1, i)
1208 rxyz(2, il) = rxyz(2, i) - alat(2)
1209 rxyz(3, il) = rxyz(3, i) - alat(3)
1210 END DO
1211
1212 END DO
1213
1214! y axis
1215 DO l2 = 0, ll2 - 1
1216
1217 in = icell(0, 0, l2, 0)
1218 icell(0, ll1, l2, ll3) = in
1219 DO ii = 1, in
1220 i = icell(ii, 0, l2, 0)
1221 il = il + 1
1222 IF (il > nn) cpabort("enlarge laymx")
1223 lay(il) = i
1224 icell(ii, ll1, l2, ll3) = il
1225 rxyz(1, il) = rxyz(1, i) + alat(1)
1226 rxyz(2, il) = rxyz(2, i)
1227 rxyz(3, il) = rxyz(3, i) + alat(3)
1228 END DO
1229
1230 in = icell(0, 0, l2, ll3 - 1)
1231 icell(0, ll1, l2, -1) = in
1232 DO ii = 1, in
1233 i = icell(ii, 0, l2, ll3 - 1)
1234 il = il + 1
1235 IF (il > nn) cpabort("enlarge laymx")
1236 lay(il) = i
1237 icell(ii, ll1, l2, -1) = il
1238 rxyz(1, il) = rxyz(1, i) + alat(1)
1239 rxyz(2, il) = rxyz(2, i)
1240 rxyz(3, il) = rxyz(3, i) - alat(3)
1241 END DO
1242
1243 in = icell(0, ll1 - 1, l2, 0)
1244 icell(0, -1, l2, ll3) = in
1245 DO ii = 1, in
1246 i = icell(ii, ll1 - 1, l2, 0)
1247 il = il + 1
1248 IF (il > nn) cpabort("enlarge laymx")
1249 lay(il) = i
1250 icell(ii, -1, l2, ll3) = il
1251 rxyz(1, il) = rxyz(1, i) - alat(1)
1252 rxyz(2, il) = rxyz(2, i)
1253 rxyz(3, il) = rxyz(3, i) + alat(3)
1254 END DO
1255
1256 in = icell(0, ll1 - 1, l2, ll3 - 1)
1257 icell(0, -1, l2, -1) = in
1258 DO ii = 1, in
1259 i = icell(ii, ll1 - 1, l2, ll3 - 1)
1260 il = il + 1
1261 IF (il > nn) cpabort("enlarge laymx")
1262 lay(il) = i
1263 icell(ii, -1, l2, -1) = il
1264 rxyz(1, il) = rxyz(1, i) - alat(1)
1265 rxyz(2, il) = rxyz(2, i)
1266 rxyz(3, il) = rxyz(3, i) - alat(3)
1267 END DO
1268
1269 END DO
1270
1271! z axis
1272 DO l3 = 0, ll3 - 1
1273
1274 in = icell(0, 0, 0, l3)
1275 icell(0, ll1, ll2, l3) = in
1276 DO ii = 1, in
1277 i = icell(ii, 0, 0, l3)
1278 il = il + 1
1279 IF (il > nn) cpabort("enlarge laymx")
1280 lay(il) = i
1281 icell(ii, ll1, ll2, l3) = il
1282 rxyz(1, il) = rxyz(1, i) + alat(1)
1283 rxyz(2, il) = rxyz(2, i) + alat(2)
1284 rxyz(3, il) = rxyz(3, i)
1285 END DO
1286
1287 in = icell(0, ll1 - 1, 0, l3)
1288 icell(0, -1, ll2, l3) = in
1289 DO ii = 1, in
1290 i = icell(ii, ll1 - 1, 0, l3)
1291 il = il + 1
1292 IF (il > nn) cpabort("enlarge laymx")
1293 lay(il) = i
1294 icell(ii, -1, ll2, l3) = il
1295 rxyz(1, il) = rxyz(1, i) - alat(1)
1296 rxyz(2, il) = rxyz(2, i) + alat(2)
1297 rxyz(3, il) = rxyz(3, i)
1298 END DO
1299
1300 in = icell(0, 0, ll2 - 1, l3)
1301 icell(0, ll1, -1, l3) = in
1302 DO ii = 1, in
1303 i = icell(ii, 0, ll2 - 1, l3)
1304 il = il + 1
1305 IF (il > nn) cpabort("enlarge laymx")
1306 lay(il) = i
1307 icell(ii, ll1, -1, l3) = il
1308 rxyz(1, il) = rxyz(1, i) + alat(1)
1309 rxyz(2, il) = rxyz(2, i) - alat(2)
1310 rxyz(3, il) = rxyz(3, i)
1311 END DO
1312
1313 in = icell(0, ll1 - 1, ll2 - 1, l3)
1314 icell(0, -1, -1, l3) = in
1315 DO ii = 1, in
1316 i = icell(ii, ll1 - 1, ll2 - 1, l3)
1317 il = il + 1
1318 IF (il > nn) cpabort("enlarge laymx")
1319 lay(il) = i
1320 icell(ii, -1, -1, l3) = il
1321 rxyz(1, il) = rxyz(1, i) - alat(1)
1322 rxyz(2, il) = rxyz(2, i) - alat(2)
1323 rxyz(3, il) = rxyz(3, i)
1324 END DO
1325
1326 END DO
1327
1328! corners
1329 in = icell(0, 0, 0, 0)
1330 icell(0, ll1, ll2, ll3) = in
1331 DO ii = 1, in
1332 i = icell(ii, 0, 0, 0)
1333 il = il + 1
1334 IF (il > nn) cpabort("enlarge laymx")
1335 lay(il) = i
1336 icell(ii, ll1, ll2, ll3) = il
1337 rxyz(1, il) = rxyz(1, i) + alat(1)
1338 rxyz(2, il) = rxyz(2, i) + alat(2)
1339 rxyz(3, il) = rxyz(3, i) + alat(3)
1340 END DO
1341
1342 in = icell(0, ll1 - 1, 0, 0)
1343 icell(0, -1, ll2, ll3) = in
1344 DO ii = 1, in
1345 i = icell(ii, ll1 - 1, 0, 0)
1346 il = il + 1
1347 IF (il > nn) cpabort("enlarge laymx")
1348 lay(il) = i
1349 icell(ii, -1, ll2, ll3) = il
1350 rxyz(1, il) = rxyz(1, i) - alat(1)
1351 rxyz(2, il) = rxyz(2, i) + alat(2)
1352 rxyz(3, il) = rxyz(3, i) + alat(3)
1353 END DO
1354
1355 in = icell(0, 0, ll2 - 1, 0)
1356 icell(0, ll1, -1, ll3) = in
1357 DO ii = 1, in
1358 i = icell(ii, 0, ll2 - 1, 0)
1359 il = il + 1
1360 IF (il > nn) cpabort("enlarge laymx")
1361 lay(il) = i
1362 icell(ii, ll1, -1, ll3) = il
1363 rxyz(1, il) = rxyz(1, i) + alat(1)
1364 rxyz(2, il) = rxyz(2, i) - alat(2)
1365 rxyz(3, il) = rxyz(3, i) + alat(3)
1366 END DO
1367
1368 in = icell(0, ll1 - 1, ll2 - 1, 0)
1369 icell(0, -1, -1, ll3) = in
1370 DO ii = 1, in
1371 i = icell(ii, ll1 - 1, ll2 - 1, 0)
1372 il = il + 1
1373 IF (il > nn) cpabort("enlarge laymx")
1374 lay(il) = i
1375 icell(ii, -1, -1, ll3) = il
1376 rxyz(1, il) = rxyz(1, i) - alat(1)
1377 rxyz(2, il) = rxyz(2, i) - alat(2)
1378 rxyz(3, il) = rxyz(3, i) + alat(3)
1379 END DO
1380
1381 in = icell(0, 0, 0, ll3 - 1)
1382 icell(0, ll1, ll2, -1) = in
1383 DO ii = 1, in
1384 i = icell(ii, 0, 0, ll3 - 1)
1385 il = il + 1
1386 IF (il > nn) cpabort("enlarge laymx")
1387 lay(il) = i
1388 icell(ii, ll1, ll2, -1) = il
1389 rxyz(1, il) = rxyz(1, i) + alat(1)
1390 rxyz(2, il) = rxyz(2, i) + alat(2)
1391 rxyz(3, il) = rxyz(3, i) - alat(3)
1392 END DO
1393
1394 in = icell(0, ll1 - 1, 0, ll3 - 1)
1395 icell(0, -1, ll2, -1) = in
1396 DO ii = 1, in
1397 i = icell(ii, ll1 - 1, 0, ll3 - 1)
1398 il = il + 1
1399 IF (il > nn) cpabort("enlarge laymx")
1400 lay(il) = i
1401 icell(ii, -1, ll2, -1) = il
1402 rxyz(1, il) = rxyz(1, i) - alat(1)
1403 rxyz(2, il) = rxyz(2, i) + alat(2)
1404 rxyz(3, il) = rxyz(3, i) - alat(3)
1405 END DO
1406
1407 in = icell(0, 0, ll2 - 1, ll3 - 1)
1408 icell(0, ll1, -1, -1) = in
1409 DO ii = 1, in
1410 i = icell(ii, 0, ll2 - 1, ll3 - 1)
1411 il = il + 1
1412 IF (il > nn) cpabort("enlarge laymx")
1413 lay(il) = i
1414 icell(ii, ll1, -1, -1) = il
1415 rxyz(1, il) = rxyz(1, i) + alat(1)
1416 rxyz(2, il) = rxyz(2, i) - alat(2)
1417 rxyz(3, il) = rxyz(3, i) - alat(3)
1418 END DO
1419
1420 in = icell(0, ll1 - 1, ll2 - 1, ll3 - 1)
1421 icell(0, -1, -1, -1) = in
1422 DO ii = 1, in
1423 i = icell(ii, ll1 - 1, ll2 - 1, ll3 - 1)
1424 il = il + 1
1425 IF (il > nn) cpabort("enlarge laymx")
1426 lay(il) = i
1427 icell(ii, -1, -1, -1) = il
1428 rxyz(1, il) = rxyz(1, i) - alat(1)
1429 rxyz(2, il) = rxyz(2, i) - alat(2)
1430 rxyz(3, il) = rxyz(3, i) - alat(3)
1431 END DO
1432
1433 ALLOCATE (lsta(2, nat))
1434 nnbrx = 12
1435 loop_nnbrx: DO
1436 ALLOCATE (lstb(nnbrx*nat), rel(5, nnbrx*nat))
1437
1438 indlstx = 0
1439
1440!$OMP PARALLEL DEFAULT(NONE) &
1441!$OMP PRIVATE(iat,cut2,iam,ii,indlst,l1,l2,l3,myspace,npr) &
1442!$OMP SHARED (indlstx,nat,nn,nnbrx,ncx,ll1,ll2,ll3,icell,lsta,lstb,lay, &
1443!$OMP rel,rxyz,cut,myspaceout)
1444
1445 npr = 1
1446!$ npr = omp_get_num_threads()
1447 iam = 0
1448!$ iam = omp_get_thread_num()
1449
1450 cut2 = cut**2
1451! assign contiguous portions of the arrays lstb and rel to the threads
1452 myspace = (nat*nnbrx)/npr
1453 IF (iam == 0) myspaceout = myspace
1454! Verlet list, relative positions
1455 indlst = 0
1456 loop_l3: DO l3 = 0, ll3 - 1
1457 loop_l2: DO l2 = 0, ll2 - 1
1458 loop_l1: DO l1 = 0, ll1 - 1
1459 loop_ii: DO ii = 1, icell(0, l1, l2, l3)
1460 iat = icell(ii, l1, l2, l3)
1461 IF (((iat - 1)*npr)/nat == iam) THEN
1462! write(*,*) 'sublstiat:iam,iat',iam,iat
1463 lsta(1, iat) = iam*myspace + indlst + 1
1464 CALL sublstiat_b(iat, nn, ncx, ll1, ll2, ll3, l1, l2, l3, myspace, &
1465 rxyz, icell, lstb(iam*myspace + 1), lay, &
1466 rel(1, iam*myspace + 1), cut2, indlst)
1467 lsta(2, iat) = iam*myspace + indlst
1468! write(*,'(a,4(x,i3),100(x,i2))') &
1469! 'iam,iat,lsta',iam,iat,lsta(1,iat),lsta(2,iat), &
1470! (lstb(j),j=lsta(1,iat),lsta(2,iat))
1471 END IF
1472 END DO loop_ii
1473 END DO loop_l1
1474 END DO loop_l2
1475 END DO loop_l3
1476!$OMP ATOMIC UPDATE
1477 indlstx = max(indlstx, indlst)
1478!$OMP END ATOMIC
1479!$OMP END PARALLEL
1480
1481 IF (indlstx >= myspaceout) THEN
1482 WRITE (10, *) count, 'NNBRX too small', nnbrx
1483 DEALLOCATE (lstb, rel)
1484 nnbrx = 3*nnbrx/2
1485 cycle loop_nnbrx
1486 END IF
1487 EXIT loop_nnbrx
1488 END DO loop_nnbrx
1489
1490 istopg = 0
1491
1492!$OMP PARALLEL DEFAULT(NONE) &
1493!$OMP PRIVATE(iam,npr,iat,iat1,iat2,lot,istop,tcoord,tcoord2, &
1494!$OMP tener,tener2,txyz,s2,s3,sz,num2,num3,numz,max_nbrs) &
1495!$OMP SHARED (nat,nnbrx,lsta,lstb,rel,ener,ener2,fxyz,coord,coord2,istopg)
1496
1497 npr = 1
1498!$ npr = omp_get_num_threads()
1499 iam = 0
1500!$ iam = omp_get_thread_num()
1501
1502 max_nbrs = 30
1503
1504 IF (npr /= 1) THEN
1505! PARALLEL CASE
1506! create temporary private scalars for reduction sum on energies and
1507! temporary private array for reduction sum on forces
1508!$OMP CRITICAL(omp_eip_bazant_silicon)
1509 ALLOCATE (txyz(3, nat), s2(max_nbrs, 8), s3(max_nbrs, 7), sz(max_nbrs, 6), &
1510 num2(max_nbrs), num3(max_nbrs), numz(max_nbrs))
1511!$OMP END CRITICAL(omp_eip_bazant_silicon)
1512 IF (iam == 0) THEN
1513 ener = 0.e0_dp
1514 ener2 = 0.e0_dp
1515 coord = 0.e0_dp
1516 coord2 = 0.e0_dp
1517 END IF
1518!$OMP DO
1519 DO iat = 1, nat
1520 fxyz(1, iat) = 0.e0_dp
1521 fxyz(2, iat) = 0.e0_dp
1522 fxyz(3, iat) = 0.e0_dp
1523 END DO
1524!$OMP BARRIER
1525
1526! Each thread treats at most lot atoms
1527 lot = int(real(nat, kind=dp)/real(npr, kind=dp) + .999999999999e0_dp)
1528 iat1 = iam*lot + 1
1529 iat2 = min((iam + 1)*lot, nat)
1530! write(*,*) 'subfeniat:iat1,iat2,iam',iat1,iat2,iam
1531 CALL subfeniat_b(iat1, iat2, nat, lsta, lstb, rel, tener, tener2, &
1532 tcoord, tcoord2, nnbrx, txyz, max_nbrs, istop, &
1533 s2(1, 1), s2(1, 2), s2(1, 3), s2(1, 4), s2(1, 5), s2(1, 6), s2(1, 7), s2(1, 8), &
1534 num2, s3(1, 1), s3(1, 2), s3(1, 3), s3(1, 4), s3(1, 5), s3(1, 6), s3(1, 7), &
1535 num3, sz(1, 1), sz(1, 2), sz(1, 3), sz(1, 4), sz(1, 5), sz(1, 6), numz)
1536
1537!$OMP CRITICAL(omp_eip_bazant_silicon)
1538 ener = ener + tener
1539 ener2 = ener2 + tener2
1540 coord = coord + tcoord
1541 coord2 = coord2 + tcoord2
1542 istopg = istopg + istop
1543 DO iat = 1, nat
1544 fxyz(1, iat) = fxyz(1, iat) + txyz(1, iat)
1545 fxyz(2, iat) = fxyz(2, iat) + txyz(2, iat)
1546 fxyz(3, iat) = fxyz(3, iat) + txyz(3, iat)
1547 END DO
1548 DEALLOCATE (txyz, s2, s3, sz, num2, num3, numz)
1549!$OMP END CRITICAL(omp_eip_bazant_silicon)
1550
1551 ELSE
1552! SERIAL CASE
1553 iat1 = 1
1554 iat2 = nat
1555 ALLOCATE (s2(max_nbrs, 8), s3(max_nbrs, 7), sz(max_nbrs, 6), &
1556 num2(max_nbrs), num3(max_nbrs), numz(max_nbrs))
1557 CALL subfeniat_b(iat1, iat2, nat, lsta, lstb, rel, ener, ener2, &
1558 coord, coord2, nnbrx, fxyz, max_nbrs, istopg, &
1559 s2(1, 1), s2(1, 2), s2(1, 3), s2(1, 4), s2(1, 5), s2(1, 6), s2(1, 7), s2(1, 8), &
1560 num2, s3(1, 1), s3(1, 2), s3(1, 3), s3(1, 4), s3(1, 5), s3(1, 6), s3(1, 7), &
1561 num3, sz(1, 1), sz(1, 2), sz(1, 3), sz(1, 4), sz(1, 5), sz(1, 6), numz)
1562 DEALLOCATE (s2, s3, sz, num2, num3, numz)
1563
1564 END IF
1565!$OMP END PARALLEL
1566
1567! write(*,*) 'ener,norm force', &
1568! ener,DNRM2(3*nat,fxyz,1)
1569 IF (istopg > 0) cpabort("DIMENSION ERROR (see WARNING above)")
1570 ener_var = ener2/nat - (ener/nat)**2
1571 coord = coord/nat
1572 coord_var = coord2/nat - coord**2
1573
1574 DEALLOCATE (rxyz, icell, lay, lsta, lstb, rel)
1575
1576 END SUBROUTINE eip_bazant_silicon
1577
1578! **************************************************************************************************
1579!> \brief ...
1580!> \param iat1 ...
1581!> \param iat2 ...
1582!> \param nat ...
1583!> \param lsta ...
1584!> \param lstb ...
1585!> \param rel ...
1586!> \param ener ...
1587!> \param ener2 ...
1588!> \param coord ...
1589!> \param coord2 ...
1590!> \param nnbrx ...
1591!> \param ff ...
1592!> \param max_nbrs ...
1593!> \param istop ...
1594!> \param s2_t0 ...
1595!> \param s2_t1 ...
1596!> \param s2_t2 ...
1597!> \param s2_t3 ...
1598!> \param s2_dx ...
1599!> \param s2_dy ...
1600!> \param s2_dz ...
1601!> \param s2_r ...
1602!> \param num2 ...
1603!> \param s3_g ...
1604!> \param s3_dg ...
1605!> \param s3_rinv ...
1606!> \param s3_dx ...
1607!> \param s3_dy ...
1608!> \param s3_dz ...
1609!> \param s3_r ...
1610!> \param num3 ...
1611!> \param sz_df ...
1612!> \param sz_sum ...
1613!> \param sz_dx ...
1614!> \param sz_dy ...
1615!> \param sz_dz ...
1616!> \param sz_r ...
1617!> \param numz ...
1618! **************************************************************************************************
1619 SUBROUTINE subfeniat_b(iat1, iat2, nat, lsta, lstb, rel, ener, ener2, &
1620 coord, coord2, nnbrx, ff, max_nbrs, istop, &
1621 s2_t0, s2_t1, s2_t2, s2_t3, s2_dx, s2_dy, s2_dz, s2_r, &
1622 num2, s3_g, s3_dg, s3_rinv, s3_dx, s3_dy, s3_dz, s3_r, &
1623 num3, sz_df, sz_sum, sz_dx, sz_dy, sz_dz, sz_r, numz)
1624! This subroutine is a modification of a subroutine that is available at
1625! http://www-math.mit.edu/~bazant/EDIP/ and for which Martin Z. Bazant
1626! and Harvard University have a 1997 copyright.
1627! The modifications were done by S. Goedecker on April 10, 2002.
1628! The routines are included with the permission of M. Bazant into this package.
1629
1630! ------------------------- VARIABLE DECLARATIONS -------------------------
1631 INTEGER :: iat1, iat2, nat, lsta(2, nat)
1632 REAL(kind=dp) :: ener, ener2, coord, coord2
1633 INTEGER :: nnbrx
1634 REAL(kind=dp) :: rel(5, nnbrx*nat)
1635 INTEGER :: lstb(nnbrx*nat)
1636 REAL(kind=dp) :: ff(3, nat)
1637 INTEGER :: max_nbrs, istop
1638 REAL(kind=dp) :: s2_t0(max_nbrs), s2_t1(max_nbrs), s2_t2(max_nbrs), s2_t3(max_nbrs), &
1639 s2_dx(max_nbrs), s2_dy(max_nbrs), s2_dz(max_nbrs), s2_r(max_nbrs)
1640 INTEGER :: num2(max_nbrs)
1641 REAL(kind=dp) :: s3_g(max_nbrs), s3_dg(max_nbrs), s3_rinv(max_nbrs), s3_dx(max_nbrs), &
1642 s3_dy(max_nbrs), s3_dz(max_nbrs), s3_r(max_nbrs)
1643 INTEGER :: num3(max_nbrs)
1644 REAL(kind=dp) :: sz_df(max_nbrs), sz_sum(max_nbrs), &
1645 sz_dx(max_nbrs), sz_dy(max_nbrs), &
1646 sz_dz(max_nbrs), sz_r(max_nbrs)
1647 INTEGER :: numz(max_nbrs)
1648
1649 INTEGER :: i, j, k, l, n, n2, n3, nj, nk, nl, nz
1650 REAL(kind=dp) :: bmc, cmbinv, coord_iat, dedrl, dedrlx, dedrly, dedrlz, den, dhdl, dhdx, &
1651 dp1, dtau, dv2dz, dv2ijx, dv2ijy, dv2ijz, dv2j, dv3dz, dv3l, dv3ljx, dv3ljy, dv3ljz, &
1652 dv3lkx, dv3lky, dv3lkz, dv3rij, dv3rijx, dv3rijy, dv3rijz, dv3rik, dv3rikx, dv3riky, &
1653 dv3rikz, dwinv, dx, dxdz, dy, dz, ener_iat, fjx, fjy, fjz, fkx, fky, fkz, fz, h, lcos, &
1654 muhalf, par_a, par_alp, par_b, par_bet, par_bg, par_c, par_cap_a, par_cap_b, par_delta, &
1655 par_eta, par_gam, par_lam, par_mu, par_palp, par_qo, par_rh, par_sig, pz, qort, r, rinv, &
1656 rmainv, rmbinv, tau, temp0, temp1, u1, u2, u3, u4, u5, winv, x, xarg
1657 REAL(kind=dp) :: xinv, xinv3, z
1658
1659! size of s2[]
1660! atom ID numbers for s2[]
1661! size of s3[]
1662! atom ID numbers for s3[]
1663! size of sz[]
1664! atom ID numbers for sz[]
1665! indices for the store arrays
1666! EDIP parameters
1667
1668 par_cap_a = 5.6714030e0_dp
1669 par_cap_b = 2.0002804e0_dp
1670 par_rh = 1.2085196e0_dp
1671 par_a = 3.1213820e0_dp
1672 par_sig = 0.5774108e0_dp
1673 par_lam = 1.4533108e0_dp
1674 par_gam = 1.1247945e0_dp
1675 par_b = 3.1213820e0_dp
1676 par_c = 2.5609104e0_dp
1677 par_delta = 78.7590539e0_dp
1678 par_mu = 0.6966326e0_dp
1679 par_qo = 312.1341346e0_dp
1680 par_palp = 1.4074424e0_dp
1681 par_bet = 0.0070975e0_dp
1682 par_alp = 3.1083847e0_dp
1683
1684 u1 = -0.165799e0_dp
1685 u2 = 32.557e0_dp
1686 u3 = 0.286198e0_dp
1687 u4 = 0.66e0_dp
1688
1689 par_bg = par_a
1690 par_eta = par_delta/par_qo
1691
1692 DO i = 1, nat
1693 ff(1, i) = 0.0e0_dp
1694 ff(2, i) = 0.0e0_dp
1695 ff(3, i) = 0.0e0_dp
1696 END DO
1697
1698 coord = 0.e0_dp
1699 coord2 = 0.e0_dp
1700 ener = 0.e0_dp
1701 ener2 = 0.e0_dp
1702 istop = 0
1703
1704! COMBINE COEFFICIENTS
1705
1706 qort = sqrt(par_qo)
1707 muhalf = par_mu*0.5e0_dp
1708 u5 = u2*u4
1709 bmc = par_b - par_c
1710 cmbinv = 1.0e0_dp/(par_c - par_b)
1711
1712! --- LEVEL 1: OUTER LOOP OVER ATOMS ---
1713
1714 atoms: DO i = iat1, iat2
1715
1716! RESET COORDINATION AND NEIGHBOR NUMBERS
1717
1718 coord_iat = 0.e0_dp
1719 ener_iat = 0.e0_dp
1720 z = 0.0e0_dp
1721 n2 = 1
1722 n3 = 1
1723 nz = 1
1724
1725! --- LEVEL 2: LOOP PREPASS OVER PAIRS ---
1726
1727 DO n = lsta(1, i), lsta(2, i)
1728 j = lstb(n)
1729
1730! PARTS OF TWO-BODY INTERACTION r<par_a
1731
1732 num2(n2) = j
1733 dx = -rel(1, n)
1734 dy = -rel(2, n)
1735 dz = -rel(3, n)
1736 r = rel(4, n)
1737 rinv = rel(5, n)
1738 rmainv = 1.e0_dp/(r - par_a)
1739 s2_t0(n2) = par_cap_a*exp(par_sig*rmainv)
1740 s2_t1(n2) = (par_cap_b*rinv)**par_rh
1741 s2_t2(n2) = par_rh*rinv
1742 s2_t3(n2) = par_sig*rmainv*rmainv
1743 s2_dx(n2) = dx
1744 s2_dy(n2) = dy
1745 s2_dz(n2) = dz
1746 s2_r(n2) = r
1747 n2 = n2 + 1
1748 IF (n2 > max_nbrs) THEN
1749 WRITE (*, *) 'WARNING enlarge max_nbrs'
1750 istop = 1
1751 RETURN
1752 END IF
1753
1754! coordination number calculated with soft cutoff between first
1755! nearest neighbor and midpoint of first and second nearest neighbor
1756 IF (r <= 2.36e0_dp) THEN
1757 coord_iat = coord_iat + 1.e0_dp
1758 ELSE IF (r >= 3.12e0_dp) THEN
1759 ELSE
1760 xarg = (r - 2.36e0_dp)*(1.e0_dp/(3.12e0_dp - 2.36e0_dp))
1761 coord_iat = coord_iat + (2*xarg + 1.e0_dp)*(xarg - 1.e0_dp)**2
1762 END IF
1763
1764! RADIAL PARTS OF THREE-BODY INTERACTION r<par_b
1765
1766 IF (r < par_bg) THEN
1767
1768 num3(n3) = j
1769 rmbinv = 1.e0_dp/(r - par_bg)
1770 temp1 = par_gam*rmbinv
1771 temp0 = exp(temp1)
1772 s3_g(n3) = temp0
1773 s3_dg(n3) = -rmbinv*temp1*temp0
1774 s3_dx(n3) = dx
1775 s3_dy(n3) = dy
1776 s3_dz(n3) = dz
1777 s3_rinv(n3) = rinv
1778 s3_r(n3) = r
1779 n3 = n3 + 1
1780 IF (n3 > max_nbrs) THEN
1781 WRITE (*, *) 'WARNING enlarge max_nbrs'
1782 istop = 1
1783 RETURN
1784 END IF
1785
1786! COORDINATION AND NEIGHBOR FUNCTION par_c<r<par_b
1787
1788 IF (r < par_b) THEN
1789 IF (r < par_c) THEN
1790 z = z + 1.e0_dp
1791 ELSE
1792 xinv = bmc/(r - par_c)
1793 xinv3 = xinv*xinv*xinv
1794 den = 1.e0_dp/(1 - xinv3)
1795 temp1 = par_alp*den
1796 fz = exp(temp1)
1797 z = z + fz
1798 numz(nz) = j
1799 sz_df(nz) = fz*temp1*den*3.e0_dp*xinv3*xinv*cmbinv
1800! df/dr
1801 sz_dx(nz) = dx
1802 sz_dy(nz) = dy
1803 sz_dz(nz) = dz
1804 sz_r(nz) = r
1805 nz = nz + 1
1806 IF (nz > max_nbrs) THEN
1807 WRITE (*, *) 'WARNING enlarge max_nbrs'
1808 istop = 1
1809 RETURN
1810 END IF
1811 END IF
1812! r < par_C
1813 END IF
1814! r < par_b
1815 END IF
1816! r < par_bg
1817 END DO
1818
1819! ZERO ACCUMULATION ARRAY FOR ENVIRONMENT FORCES
1820
1821 DO nl = 1, nz - 1
1822 sz_sum(nl) = 0.e0_dp
1823 END DO
1824
1825! ENVIRONMENT-DEPENDENCE OF PAIR INTERACTION
1826
1827 temp0 = par_bet*z
1828 pz = par_palp*exp(-temp0*z)
1829! bond order
1830 dp1 = -2.e0_dp*temp0*pz
1831! derivative of bond order
1832
1833! --- LEVEL 2: LOOP FOR PAIR INTERACTIONS ---
1834
1835 DO nj = 1, n2 - 1
1836
1837 temp0 = s2_t1(nj) - pz
1838
1839! two-body energy V2(rij,Z)
1840
1841 ener_iat = ener_iat + temp0*s2_t0(nj)
1842
1843! two-body forces
1844
1845 dv2j = -s2_t0(nj)*(s2_t1(nj)*s2_t2(nj) + temp0*s2_t3(nj))
1846! dV2/dr
1847 dv2ijx = dv2j*s2_dx(nj)
1848 dv2ijy = dv2j*s2_dy(nj)
1849 dv2ijz = dv2j*s2_dz(nj)
1850 ff(1, i) = ff(1, i) + dv2ijx
1851 ff(2, i) = ff(2, i) + dv2ijy
1852 ff(3, i) = ff(3, i) + dv2ijz
1853 j = num2(nj)
1854 ff(1, j) = ff(1, j) - dv2ijx
1855 ff(2, j) = ff(2, j) - dv2ijy
1856 ff(3, j) = ff(3, j) - dv2ijz
1857
1858! --- LEVEL 3: LOOP FOR PAIR COORDINATION FORCES ---
1859
1860 dv2dz = -dp1*s2_t0(nj)
1861 DO nl = 1, nz - 1
1862 sz_sum(nl) = sz_sum(nl) + dv2dz
1863 END DO
1864
1865 END DO
1866
1867! COORDINATION-DEPENDENCE OF THREE-BODY INTERACTION
1868
1869 winv = qort*exp(-muhalf*z)
1870! inverse width of angular function
1871 dwinv = -muhalf*winv
1872! its derivative
1873 temp0 = exp(-u4*z)
1874 tau = u1 + u2*temp0*(u3 - temp0)
1875! -cosine of angular minimum
1876 dtau = u5*temp0*(2*temp0 - u3)
1877! its derivative
1878
1879! --- LEVEL 2: FIRST LOOP FOR THREE-BODY INTERACTIONS ---
1880
1881 DO nj = 1, n3 - 2
1882
1883 j = num3(nj)
1884
1885! --- LEVEL 3: SECOND LOOP FOR THREE-BODY INTERACTIONS ---
1886
1887 DO nk = nj + 1, n3 - 1
1888
1889 k = num3(nk)
1890
1891! angular function h(l,Z)
1892
1893 lcos = s3_dx(nj)*s3_dx(nk) + s3_dy(nj)*s3_dy(nk) + s3_dz(nj)*s3_dz(nk)
1894 x = (lcos + tau)*winv
1895 temp0 = exp(-x*x)
1896
1897 h = par_lam*(1 - temp0 + par_eta*x*x)
1898 dhdx = 2*par_lam*x*(temp0 + par_eta)
1899
1900 dhdl = dhdx*winv
1901
1902! three-body energy
1903
1904 temp1 = s3_g(nj)*s3_g(nk)
1905 ener_iat = ener_iat + temp1*h
1906
1907! (-) radial force on atom j
1908
1909 dv3rij = s3_dg(nj)*s3_g(nk)*h
1910 dv3rijx = dv3rij*s3_dx(nj)
1911 dv3rijy = dv3rij*s3_dy(nj)
1912 dv3rijz = dv3rij*s3_dz(nj)
1913 fjx = dv3rijx
1914 fjy = dv3rijy
1915 fjz = dv3rijz
1916
1917! (-) radial force on atom k
1918
1919 dv3rik = s3_g(nj)*s3_dg(nk)*h
1920 dv3rikx = dv3rik*s3_dx(nk)
1921 dv3riky = dv3rik*s3_dy(nk)
1922 dv3rikz = dv3rik*s3_dz(nk)
1923 fkx = dv3rikx
1924 fky = dv3riky
1925 fkz = dv3rikz
1926
1927! (-) angular force on j
1928
1929 dv3l = temp1*dhdl
1930 dv3ljx = dv3l*(s3_dx(nk) - lcos*s3_dx(nj))*s3_rinv(nj)
1931 dv3ljy = dv3l*(s3_dy(nk) - lcos*s3_dy(nj))*s3_rinv(nj)
1932 dv3ljz = dv3l*(s3_dz(nk) - lcos*s3_dz(nj))*s3_rinv(nj)
1933 fjx = fjx + dv3ljx
1934 fjy = fjy + dv3ljy
1935 fjz = fjz + dv3ljz
1936
1937! (-) angular force on k
1938
1939 dv3lkx = dv3l*(s3_dx(nj) - lcos*s3_dx(nk))*s3_rinv(nk)
1940 dv3lky = dv3l*(s3_dy(nj) - lcos*s3_dy(nk))*s3_rinv(nk)
1941 dv3lkz = dv3l*(s3_dz(nj) - lcos*s3_dz(nk))*s3_rinv(nk)
1942 fkx = fkx + dv3lkx
1943 fky = fky + dv3lky
1944 fkz = fkz + dv3lkz
1945
1946! apply radial + angular forces to i, j, k
1947
1948 ff(1, j) = ff(1, j) - fjx
1949 ff(2, j) = ff(2, j) - fjy
1950 ff(3, j) = ff(3, j) - fjz
1951 ff(1, k) = ff(1, k) - fkx
1952 ff(2, k) = ff(2, k) - fky
1953 ff(3, k) = ff(3, k) - fkz
1954 ff(1, i) = ff(1, i) + fjx + fkx
1955 ff(2, i) = ff(2, i) + fjy + fky
1956 ff(3, i) = ff(3, i) + fjz + fkz
1957
1958! prefactor for 4-body forces from coordination
1959 dxdz = dwinv*(lcos + tau) + winv*dtau
1960 dv3dz = temp1*dhdx*dxdz
1961
1962! --- LEVEL 4: LOOP FOR THREE-BODY COORDINATION FORCES ---
1963
1964 DO nl = 1, nz - 1
1965 sz_sum(nl) = sz_sum(nl) + dv3dz
1966 END DO
1967 END DO
1968 END DO
1969
1970! --- LEVEL 2: LOOP TO APPLY COORDINATION FORCES ---
1971
1972 DO nl = 1, nz - 1
1973
1974 dedrl = sz_sum(nl)*sz_df(nl)
1975 dedrlx = dedrl*sz_dx(nl)
1976 dedrly = dedrl*sz_dy(nl)
1977 dedrlz = dedrl*sz_dz(nl)
1978 ff(1, i) = ff(1, i) + dedrlx
1979 ff(2, i) = ff(2, i) + dedrly
1980 ff(3, i) = ff(3, i) + dedrlz
1981 l = numz(nl)
1982 ff(1, l) = ff(1, l) - dedrlx
1983 ff(2, l) = ff(2, l) - dedrly
1984 ff(3, l) = ff(3, l) - dedrlz
1985
1986 END DO
1987
1988 coord = coord + coord_iat
1989 coord2 = coord2 + coord_iat**2
1990 ener = ener + ener_iat
1991 ener2 = ener2 + ener_iat**2
1992
1993 END DO atoms
1994
1995 RETURN
1996 END SUBROUTINE subfeniat_b
1997
1998! **************************************************************************************************
1999!> \brief ...
2000!> \param iat ...
2001!> \param nn ...
2002!> \param ncx ...
2003!> \param ll1 ...
2004!> \param ll2 ...
2005!> \param ll3 ...
2006!> \param l1 ...
2007!> \param l2 ...
2008!> \param l3 ...
2009!> \param myspace ...
2010!> \param rxyz ...
2011!> \param icell ...
2012!> \param lstb ...
2013!> \param lay ...
2014!> \param rel ...
2015!> \param cut2 ...
2016!> \param indlst ...
2017! **************************************************************************************************
2018 SUBROUTINE sublstiat_b(iat, nn, ncx, ll1, ll2, ll3, l1, l2, l3, myspace, &
2019 rxyz, icell, lstb, lay, rel, cut2, indlst)
2020! finds the neighbours of atom iat (specified by lsta and lstb) and and
2021! the relative position rel of iat with respect to these neighbours
2022 INTEGER :: iat, nn, ncx, ll1, ll2, ll3, l1, l2, l3, &
2023 myspace
2024 REAL(kind=dp) :: rxyz
2025 INTEGER :: icell, lstb, lay
2026 REAL(kind=dp) :: rel, cut2
2027 INTEGER :: indlst
2028
2029 dimension rxyz(3, nn), lay(nn), icell(0:ncx, -1:ll1, -1:ll2, -1:ll3), &
2030 lstb(0:myspace - 1), rel(5, 0:myspace - 1)
2031
2032 INTEGER :: jat, k1, k2, k3, jj
2033 REAL(kind=dp) :: rr2, tt, tti, xrel, yrel, zrel
2034
2035 DO k3 = l3 - 1, l3 + 1
2036 DO k2 = l2 - 1, l2 + 1
2037 DO k1 = l1 - 1, l1 + 1
2038 DO jj = 1, icell(0, k1, k2, k3)
2039 jat = icell(jj, k1, k2, k3)
2040 IF (jat == iat) cycle
2041 xrel = rxyz(1, iat) - rxyz(1, jat)
2042 yrel = rxyz(2, iat) - rxyz(2, jat)
2043 zrel = rxyz(3, iat) - rxyz(3, jat)
2044 rr2 = xrel**2 + yrel**2 + zrel**2
2045 IF (rr2 <= cut2) THEN
2046 indlst = min(indlst, myspace - 1)
2047 lstb(indlst) = lay(jat)
2048! write(*,*) 'iat,indlst,lay(jat)',iat,indlst,lay(jat)
2049 tt = sqrt(rr2)
2050 tti = 1.e0_dp/tt
2051 rel(1, indlst) = xrel*tti
2052 rel(2, indlst) = yrel*tti
2053 rel(3, indlst) = zrel*tti
2054 rel(4, indlst) = tt
2055 rel(5, indlst) = tti
2056 indlst = indlst + 1
2057 END IF
2058 END DO
2059 END DO
2060 END DO
2061 END DO
2062
2063 RETURN
2064 END SUBROUTINE sublstiat_b
2065
2066! **************************************************************************************************
2067!> \brief Lenosky's "highly optimized empirical potential model of silicon"
2068!> by Stefan Goedecker
2069!> \param nat number of atoms
2070!> \param alat lattice constants of the orthorombic box containing the particles
2071!> \param rxyz0 atomic positions in Angstrom, may be modified on output.
2072!> If an atom is outside the box the program will bring it back
2073!> into the box by translations through alat
2074!> \param fxyz forces in eV/A
2075!> \param ener total energy in eV
2076!> \param coord average coordination number
2077!> \param ener_var variance of the energy/atom
2078!> \param coord_var variance of the coordination number
2079!> \param count count is increased by one per call, has to be initialized
2080!> to 0.e0_dp before first call of eip_bazant
2081!> \par Literature
2082!> T. Lenosky, et. al.: Highly optimized empirical potential model of silicon;
2083!> Modeling Simul. Sci. Eng., 8 (2000)
2084!> S. Goedecker: Optimization and parallelization of a force field for silicon
2085!> using OpenMP; CPC 148, 1 (2002)
2086!> \par History
2087!> 03.2006 initial create [tdk]
2088!> \author Thomas D. Kuehne (tkuehne@cp2k.org)
2089! **************************************************************************************************
2090 SUBROUTINE eip_lenosky_silicon(nat, alat, rxyz0, fxyz, ener, coord, ener_var, &
2091 coord_var, count)
2092
2093 INTEGER :: nat
2094 REAL(kind=dp) :: alat, rxyz0, fxyz, ener, coord, &
2095 ener_var, coord_var, count
2096
2097 dimension rxyz0(3, nat), fxyz(3, nat), alat(3)
2098 REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :) :: rxyz
2099 INTEGER, ALLOCATABLE, DIMENSION(:, :) :: lsta
2100 INTEGER, ALLOCATABLE, DIMENSION(:) :: lstb
2101 INTEGER, ALLOCATABLE, DIMENSION(:) :: lay
2102 INTEGER, ALLOCATABLE, DIMENSION(:, :, :, :) :: icell
2103 REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :) :: rel
2104 REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :) :: txyz
2105 REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :) :: f2ij, f3ij, f3ik
2106
2107 REAL(kind=dp) :: coord2, cut, cut2, ener2, tcoord, &
2108 tcoord2, tener, tener2
2109 INTEGER :: i, iam, iat, iat1, iat2, ii, il, in, indlst, indlstx, &
2110 istop, istopg, l1, l2, l3, ll1, ll2, ll3, lot, ncx, nn, &
2111 nnbrx, npjkx, npjx, laymx, npr, rlc1i, rlc2i, rlc3i, &
2112 myspace, myspaceout
2113
2114! tmax_phi= 0.4500000e+01_dp
2115! cut=tmax_phi
2116 cut = 0.4500000e+01_dp
2117
2118 IF (count == 0) OPEN (unit=10, file='lenosky.mon', status='unknown')
2119 count = count + 1.e0_dp
2120
2121! linear scaling calculation of verlet list
2122 ll1 = int(alat(1)/cut)
2123 IF (ll1 < 1) cpabort("alat(1) too small")
2124 ll2 = int(alat(2)/cut)
2125 IF (ll2 < 1) cpabort("alat(2) too small")
2126 ll3 = int(alat(3)/cut)
2127 IF (ll3 < 1) cpabort("alat(3) too small")
2128
2129! determine number of threads
2130 npr = 1
2131!$OMP PARALLEL PRIVATE(iam) SHARED (npr) DEFAULT(NONE)
2132!$ iam = omp_get_thread_num()
2133!$ if (iam .eq. 0) npr = omp_get_num_threads()
2134!$OMP END PARALLEL
2135
2136! linear scaling calculation of verlet list
2137
2138 IF (npr <= 1) THEN !serial if too few processors to gain by parallelizing
2139
2140! set ncx for serial case, ncx for parallel case set below
2141 ncx = 16
2142 loop_ncx_s: DO
2143 ALLOCATE (icell(0:ncx, -1:ll1, -1:ll2, -1:ll3))
2144 icell(0, -1:ll1, -1:ll2, -1:ll3) = 0
2145 rlc1i = int(ll1/alat(1))
2146 rlc2i = int(ll2/alat(2))
2147 rlc3i = int(ll3/alat(3))
2148
2149 loop_iat_s: DO iat = 1, nat
2150 rxyz0(1, iat) = modulo(modulo(rxyz0(1, iat), alat(1)), alat(1))
2151 rxyz0(2, iat) = modulo(modulo(rxyz0(2, iat), alat(2)), alat(2))
2152 rxyz0(3, iat) = modulo(modulo(rxyz0(3, iat), alat(3)), alat(3))
2153 l1 = int(rxyz0(1, iat)*rlc1i)
2154 l2 = int(rxyz0(2, iat)*rlc2i)
2155 l3 = int(rxyz0(3, iat)*rlc3i)
2156
2157 ii = icell(0, l1, l2, l3)
2158 ii = ii + 1
2159 icell(0, l1, l2, l3) = ii
2160 IF (ii > ncx) THEN
2161 WRITE (10, *) count, 'NCX too small', ncx
2162 DEALLOCATE (icell)
2163 ncx = ncx*2
2164 cycle loop_ncx_s
2165 END IF
2166 icell(ii, l1, l2, l3) = iat
2167 END DO loop_iat_s
2168 EXIT loop_ncx_s
2169 END DO loop_ncx_s
2170
2171 ELSE ! parallel case
2172
2173! periodization of particles can be done in parallel
2174!$OMP PARALLEL DO SHARED (alat,nat,rxyz0) PRIVATE(iat) DEFAULT(NONE)
2175 DO iat = 1, nat
2176 rxyz0(1, iat) = modulo(modulo(rxyz0(1, iat), alat(1)), alat(1))
2177 rxyz0(2, iat) = modulo(modulo(rxyz0(2, iat), alat(2)), alat(2))
2178 rxyz0(3, iat) = modulo(modulo(rxyz0(3, iat), alat(3)), alat(3))
2179 END DO
2180!$OMP END PARALLEL DO
2181
2182! assignment to cell is done serially
2183! set ncx for parallel case, ncx for serial case set above
2184 ncx = 16
2185 loop_ncx_p: DO
2186 ALLOCATE (icell(0:ncx, -1:ll1, -1:ll2, -1:ll3))
2187 icell(0, -1:ll1, -1:ll2, -1:ll3) = 0
2188 rlc1i = int(ll1/alat(1))
2189 rlc2i = int(ll2/alat(2))
2190 rlc3i = int(ll3/alat(3))
2191
2192 loop_iat_p: DO iat = 1, nat
2193 l1 = int(rxyz0(1, iat)*rlc1i)
2194 l2 = int(rxyz0(2, iat)*rlc2i)
2195 l3 = int(rxyz0(3, iat)*rlc3i)
2196 ii = icell(0, l1, l2, l3)
2197 ii = ii + 1
2198 icell(0, l1, l2, l3) = ii
2199 IF (ii > ncx) THEN
2200 WRITE (10, *) count, 'NCX too small', ncx
2201 DEALLOCATE (icell)
2202 ncx = ncx*2
2203 cycle loop_ncx_p
2204 END IF
2205 icell(ii, l1, l2, l3) = iat
2206 END DO loop_iat_p
2207 EXIT loop_ncx_p
2208 END DO loop_ncx_p
2209
2210 END IF
2211
2212! duplicate all atoms within boundary layer
2213 laymx = ncx*(2*ll1*ll2 + 2*ll1*ll3 + 2*ll2*ll3 + 4*ll1 + 4*ll2 + 4*ll3 + 8)
2214 nn = nat + laymx
2215 ALLOCATE (rxyz(3, nn), lay(nn))
2216 DO iat = 1, nat
2217 lay(iat) = iat
2218 rxyz(1, iat) = rxyz0(1, iat)
2219 rxyz(2, iat) = rxyz0(2, iat)
2220 rxyz(3, iat) = rxyz0(3, iat)
2221 END DO
2222 il = nat
2223! xy plane
2224 DO l2 = 0, ll2 - 1
2225 DO l1 = 0, ll1 - 1
2226
2227 in = icell(0, l1, l2, 0)
2228 icell(0, l1, l2, ll3) = in
2229 DO ii = 1, in
2230 i = icell(ii, l1, l2, 0)
2231 il = il + 1
2232 IF (il > nn) cpabort("enlarge laymx")
2233 lay(il) = i
2234 icell(ii, l1, l2, ll3) = il
2235 rxyz(1, il) = rxyz(1, i)
2236 rxyz(2, il) = rxyz(2, i)
2237 rxyz(3, il) = rxyz(3, i) + alat(3)
2238 END DO
2239
2240 in = icell(0, l1, l2, ll3 - 1)
2241 icell(0, l1, l2, -1) = in
2242 DO ii = 1, in
2243 i = icell(ii, l1, l2, ll3 - 1)
2244 il = il + 1
2245 IF (il > nn) cpabort("enlarge laymx")
2246 lay(il) = i
2247 icell(ii, l1, l2, -1) = il
2248 rxyz(1, il) = rxyz(1, i)
2249 rxyz(2, il) = rxyz(2, i)
2250 rxyz(3, il) = rxyz(3, i) - alat(3)
2251 END DO
2252
2253 END DO
2254 END DO
2255
2256! yz plane
2257 DO l3 = 0, ll3 - 1
2258 DO l2 = 0, ll2 - 1
2259
2260 in = icell(0, 0, l2, l3)
2261 icell(0, ll1, l2, l3) = in
2262 DO ii = 1, in
2263 i = icell(ii, 0, l2, l3)
2264 il = il + 1
2265 IF (il > nn) cpabort("enlarge laymx")
2266 lay(il) = i
2267 icell(ii, ll1, l2, l3) = il
2268 rxyz(1, il) = rxyz(1, i) + alat(1)
2269 rxyz(2, il) = rxyz(2, i)
2270 rxyz(3, il) = rxyz(3, i)
2271 END DO
2272
2273 in = icell(0, ll1 - 1, l2, l3)
2274 icell(0, -1, l2, l3) = in
2275 DO ii = 1, in
2276 i = icell(ii, ll1 - 1, l2, l3)
2277 il = il + 1
2278 IF (il > nn) cpabort("enlarge laymx")
2279 lay(il) = i
2280 icell(ii, -1, l2, l3) = il
2281 rxyz(1, il) = rxyz(1, i) - alat(1)
2282 rxyz(2, il) = rxyz(2, i)
2283 rxyz(3, il) = rxyz(3, i)
2284 END DO
2285
2286 END DO
2287 END DO
2288
2289! xz plane
2290 DO l3 = 0, ll3 - 1
2291 DO l1 = 0, ll1 - 1
2292
2293 in = icell(0, l1, 0, l3)
2294 icell(0, l1, ll2, l3) = in
2295 DO ii = 1, in
2296 i = icell(ii, l1, 0, l3)
2297 il = il + 1
2298 IF (il > nn) cpabort("enlarge laymx")
2299 lay(il) = i
2300 icell(ii, l1, ll2, l3) = il
2301 rxyz(1, il) = rxyz(1, i)
2302 rxyz(2, il) = rxyz(2, i) + alat(2)
2303 rxyz(3, il) = rxyz(3, i)
2304 END DO
2305
2306 in = icell(0, l1, ll2 - 1, l3)
2307 icell(0, l1, -1, l3) = in
2308 DO ii = 1, in
2309 i = icell(ii, l1, ll2 - 1, l3)
2310 il = il + 1
2311 IF (il > nn) cpabort("enlarge laymx")
2312 lay(il) = i
2313 icell(ii, l1, -1, l3) = il
2314 rxyz(1, il) = rxyz(1, i)
2315 rxyz(2, il) = rxyz(2, i) - alat(2)
2316 rxyz(3, il) = rxyz(3, i)
2317 END DO
2318
2319 END DO
2320 END DO
2321
2322! x axis
2323 DO l1 = 0, ll1 - 1
2324
2325 in = icell(0, l1, 0, 0)
2326 icell(0, l1, ll2, ll3) = in
2327 DO ii = 1, in
2328 i = icell(ii, l1, 0, 0)
2329 il = il + 1
2330 IF (il > nn) cpabort("enlarge laymx")
2331 lay(il) = i
2332 icell(ii, l1, ll2, ll3) = il
2333 rxyz(1, il) = rxyz(1, i)
2334 rxyz(2, il) = rxyz(2, i) + alat(2)
2335 rxyz(3, il) = rxyz(3, i) + alat(3)
2336 END DO
2337
2338 in = icell(0, l1, 0, ll3 - 1)
2339 icell(0, l1, ll2, -1) = in
2340 DO ii = 1, in
2341 i = icell(ii, l1, 0, ll3 - 1)
2342 il = il + 1
2343 IF (il > nn) cpabort("enlarge laymx")
2344 lay(il) = i
2345 icell(ii, l1, ll2, -1) = il
2346 rxyz(1, il) = rxyz(1, i)
2347 rxyz(2, il) = rxyz(2, i) + alat(2)
2348 rxyz(3, il) = rxyz(3, i) - alat(3)
2349 END DO
2350
2351 in = icell(0, l1, ll2 - 1, 0)
2352 icell(0, l1, -1, ll3) = in
2353 DO ii = 1, in
2354 i = icell(ii, l1, ll2 - 1, 0)
2355 il = il + 1
2356 IF (il > nn) cpabort("enlarge laymx")
2357 lay(il) = i
2358 icell(ii, l1, -1, ll3) = il
2359 rxyz(1, il) = rxyz(1, i)
2360 rxyz(2, il) = rxyz(2, i) - alat(2)
2361 rxyz(3, il) = rxyz(3, i) + alat(3)
2362 END DO
2363
2364 in = icell(0, l1, ll2 - 1, ll3 - 1)
2365 icell(0, l1, -1, -1) = in
2366 DO ii = 1, in
2367 i = icell(ii, l1, ll2 - 1, ll3 - 1)
2368 il = il + 1
2369 IF (il > nn) cpabort("enlarge laymx")
2370 lay(il) = i
2371 icell(ii, l1, -1, -1) = il
2372 rxyz(1, il) = rxyz(1, i)
2373 rxyz(2, il) = rxyz(2, i) - alat(2)
2374 rxyz(3, il) = rxyz(3, i) - alat(3)
2375 END DO
2376
2377 END DO
2378
2379! y axis
2380 DO l2 = 0, ll2 - 1
2381
2382 in = icell(0, 0, l2, 0)
2383 icell(0, ll1, l2, ll3) = in
2384 DO ii = 1, in
2385 i = icell(ii, 0, l2, 0)
2386 il = il + 1
2387 IF (il > nn) cpabort("enlarge laymx")
2388 lay(il) = i
2389 icell(ii, ll1, l2, ll3) = il
2390 rxyz(1, il) = rxyz(1, i) + alat(1)
2391 rxyz(2, il) = rxyz(2, i)
2392 rxyz(3, il) = rxyz(3, i) + alat(3)
2393 END DO
2394
2395 in = icell(0, 0, l2, ll3 - 1)
2396 icell(0, ll1, l2, -1) = in
2397 DO ii = 1, in
2398 i = icell(ii, 0, l2, ll3 - 1)
2399 il = il + 1
2400 IF (il > nn) cpabort("enlarge laymx")
2401 lay(il) = i
2402 icell(ii, ll1, l2, -1) = il
2403 rxyz(1, il) = rxyz(1, i) + alat(1)
2404 rxyz(2, il) = rxyz(2, i)
2405 rxyz(3, il) = rxyz(3, i) - alat(3)
2406 END DO
2407
2408 in = icell(0, ll1 - 1, l2, 0)
2409 icell(0, -1, l2, ll3) = in
2410 DO ii = 1, in
2411 i = icell(ii, ll1 - 1, l2, 0)
2412 il = il + 1
2413 IF (il > nn) cpabort("enlarge laymx")
2414 lay(il) = i
2415 icell(ii, -1, l2, ll3) = il
2416 rxyz(1, il) = rxyz(1, i) - alat(1)
2417 rxyz(2, il) = rxyz(2, i)
2418 rxyz(3, il) = rxyz(3, i) + alat(3)
2419 END DO
2420
2421 in = icell(0, ll1 - 1, l2, ll3 - 1)
2422 icell(0, -1, l2, -1) = in
2423 DO ii = 1, in
2424 i = icell(ii, ll1 - 1, l2, ll3 - 1)
2425 il = il + 1
2426 IF (il > nn) cpabort("enlarge laymx")
2427 lay(il) = i
2428 icell(ii, -1, l2, -1) = il
2429 rxyz(1, il) = rxyz(1, i) - alat(1)
2430 rxyz(2, il) = rxyz(2, i)
2431 rxyz(3, il) = rxyz(3, i) - alat(3)
2432 END DO
2433
2434 END DO
2435
2436! z axis
2437 DO l3 = 0, ll3 - 1
2438
2439 in = icell(0, 0, 0, l3)
2440 icell(0, ll1, ll2, l3) = in
2441 DO ii = 1, in
2442 i = icell(ii, 0, 0, l3)
2443 il = il + 1
2444 IF (il > nn) cpabort("enlarge laymx")
2445 lay(il) = i
2446 icell(ii, ll1, ll2, l3) = il
2447 rxyz(1, il) = rxyz(1, i) + alat(1)
2448 rxyz(2, il) = rxyz(2, i) + alat(2)
2449 rxyz(3, il) = rxyz(3, i)
2450 END DO
2451
2452 in = icell(0, ll1 - 1, 0, l3)
2453 icell(0, -1, ll2, l3) = in
2454 DO ii = 1, in
2455 i = icell(ii, ll1 - 1, 0, l3)
2456 il = il + 1
2457 IF (il > nn) cpabort("enlarge laymx")
2458 lay(il) = i
2459 icell(ii, -1, ll2, l3) = il
2460 rxyz(1, il) = rxyz(1, i) - alat(1)
2461 rxyz(2, il) = rxyz(2, i) + alat(2)
2462 rxyz(3, il) = rxyz(3, i)
2463 END DO
2464
2465 in = icell(0, 0, ll2 - 1, l3)
2466 icell(0, ll1, -1, l3) = in
2467 DO ii = 1, in
2468 i = icell(ii, 0, ll2 - 1, l3)
2469 il = il + 1
2470 IF (il > nn) cpabort("enlarge laymx")
2471 lay(il) = i
2472 icell(ii, ll1, -1, l3) = il
2473 rxyz(1, il) = rxyz(1, i) + alat(1)
2474 rxyz(2, il) = rxyz(2, i) - alat(2)
2475 rxyz(3, il) = rxyz(3, i)
2476 END DO
2477
2478 in = icell(0, ll1 - 1, ll2 - 1, l3)
2479 icell(0, -1, -1, l3) = in
2480 DO ii = 1, in
2481 i = icell(ii, ll1 - 1, ll2 - 1, l3)
2482 il = il + 1
2483 IF (il > nn) cpabort("enlarge laymx")
2484 lay(il) = i
2485 icell(ii, -1, -1, l3) = il
2486 rxyz(1, il) = rxyz(1, i) - alat(1)
2487 rxyz(2, il) = rxyz(2, i) - alat(2)
2488 rxyz(3, il) = rxyz(3, i)
2489 END DO
2490
2491 END DO
2492
2493! corners
2494 in = icell(0, 0, 0, 0)
2495 icell(0, ll1, ll2, ll3) = in
2496 DO ii = 1, in
2497 i = icell(ii, 0, 0, 0)
2498 il = il + 1
2499 IF (il > nn) cpabort("enlarge laymx")
2500 lay(il) = i
2501 icell(ii, ll1, ll2, ll3) = il
2502 rxyz(1, il) = rxyz(1, i) + alat(1)
2503 rxyz(2, il) = rxyz(2, i) + alat(2)
2504 rxyz(3, il) = rxyz(3, i) + alat(3)
2505 END DO
2506
2507 in = icell(0, ll1 - 1, 0, 0)
2508 icell(0, -1, ll2, ll3) = in
2509 DO ii = 1, in
2510 i = icell(ii, ll1 - 1, 0, 0)
2511 il = il + 1
2512 IF (il > nn) cpabort("enlarge laymx")
2513 lay(il) = i
2514 icell(ii, -1, ll2, ll3) = il
2515 rxyz(1, il) = rxyz(1, i) - alat(1)
2516 rxyz(2, il) = rxyz(2, i) + alat(2)
2517 rxyz(3, il) = rxyz(3, i) + alat(3)
2518 END DO
2519
2520 in = icell(0, 0, ll2 - 1, 0)
2521 icell(0, ll1, -1, ll3) = in
2522 DO ii = 1, in
2523 i = icell(ii, 0, ll2 - 1, 0)
2524 il = il + 1
2525 IF (il > nn) cpabort("enlarge laymx")
2526 lay(il) = i
2527 icell(ii, ll1, -1, ll3) = il
2528 rxyz(1, il) = rxyz(1, i) + alat(1)
2529 rxyz(2, il) = rxyz(2, i) - alat(2)
2530 rxyz(3, il) = rxyz(3, i) + alat(3)
2531 END DO
2532
2533 in = icell(0, ll1 - 1, ll2 - 1, 0)
2534 icell(0, -1, -1, ll3) = in
2535 DO ii = 1, in
2536 i = icell(ii, ll1 - 1, ll2 - 1, 0)
2537 il = il + 1
2538 IF (il > nn) cpabort("enlarge laymx")
2539 lay(il) = i
2540 icell(ii, -1, -1, ll3) = il
2541 rxyz(1, il) = rxyz(1, i) - alat(1)
2542 rxyz(2, il) = rxyz(2, i) - alat(2)
2543 rxyz(3, il) = rxyz(3, i) + alat(3)
2544 END DO
2545
2546 in = icell(0, 0, 0, ll3 - 1)
2547 icell(0, ll1, ll2, -1) = in
2548 DO ii = 1, in
2549 i = icell(ii, 0, 0, ll3 - 1)
2550 il = il + 1
2551 IF (il > nn) cpabort("enlarge laymx")
2552 lay(il) = i
2553 icell(ii, ll1, ll2, -1) = il
2554 rxyz(1, il) = rxyz(1, i) + alat(1)
2555 rxyz(2, il) = rxyz(2, i) + alat(2)
2556 rxyz(3, il) = rxyz(3, i) - alat(3)
2557 END DO
2558
2559 in = icell(0, ll1 - 1, 0, ll3 - 1)
2560 icell(0, -1, ll2, -1) = in
2561 DO ii = 1, in
2562 i = icell(ii, ll1 - 1, 0, ll3 - 1)
2563 il = il + 1
2564 IF (il > nn) cpabort("enlarge laymx")
2565 lay(il) = i
2566 icell(ii, -1, ll2, -1) = il
2567 rxyz(1, il) = rxyz(1, i) - alat(1)
2568 rxyz(2, il) = rxyz(2, i) + alat(2)
2569 rxyz(3, il) = rxyz(3, i) - alat(3)
2570 END DO
2571
2572 in = icell(0, 0, ll2 - 1, ll3 - 1)
2573 icell(0, ll1, -1, -1) = in
2574 DO ii = 1, in
2575 i = icell(ii, 0, ll2 - 1, ll3 - 1)
2576 il = il + 1
2577 IF (il > nn) cpabort("enlarge laymx")
2578 lay(il) = i
2579 icell(ii, ll1, -1, -1) = il
2580 rxyz(1, il) = rxyz(1, i) + alat(1)
2581 rxyz(2, il) = rxyz(2, i) - alat(2)
2582 rxyz(3, il) = rxyz(3, i) - alat(3)
2583 END DO
2584
2585 in = icell(0, ll1 - 1, ll2 - 1, ll3 - 1)
2586 icell(0, -1, -1, -1) = in
2587 DO ii = 1, in
2588 i = icell(ii, ll1 - 1, ll2 - 1, ll3 - 1)
2589 il = il + 1
2590 IF (il > nn) cpabort("enlarge laymx")
2591 lay(il) = i
2592 icell(ii, -1, -1, -1) = il
2593 rxyz(1, il) = rxyz(1, i) - alat(1)
2594 rxyz(2, il) = rxyz(2, i) - alat(2)
2595 rxyz(3, il) = rxyz(3, i) - alat(3)
2596 END DO
2597
2598 ALLOCATE (lsta(2, nat))
2599 nnbrx = 36
2600 loop_nnbrx: DO
2601 ALLOCATE (lstb(nnbrx*nat), rel(5, nnbrx*nat))
2602
2603 indlstx = 0
2604
2605!$OMP PARALLEL DEFAULT(NONE) &
2606!$OMP PRIVATE(iat,cut2,iam,ii,indlst,l1,l2,l3,myspace,npr) &
2607!$OMP SHARED (indlstx,nat,nn,nnbrx,ncx,ll1,ll2,ll3,icell,lsta,lstb,lay, &
2608!$OMP rel,rxyz,cut,myspaceout)
2609
2610 npr = 1
2611!$ npr = omp_get_num_threads()
2612 iam = 0
2613!$ iam = omp_get_thread_num()
2614
2615 cut2 = cut**2
2616! assign contiguous portions of the arrays lstb and rel to the threads
2617 myspace = (nat*nnbrx)/npr
2618 IF (iam == 0) myspaceout = myspace
2619! Verlet list, relative positions
2620 indlst = 0
2621 loop_l3: DO l3 = 0, ll3 - 1
2622 loop_l2: DO l2 = 0, ll2 - 1
2623 loop_l1: DO l1 = 0, ll1 - 1
2624 loop_ii: DO ii = 1, icell(0, l1, l2, l3)
2625 iat = icell(ii, l1, l2, l3)
2626 IF (((iat - 1)*npr)/nat == iam) THEN
2627! write(*,*) 'sublstiat:iam,iat',iam,iat
2628 lsta(1, iat) = iam*myspace + indlst + 1
2629 CALL sublstiat_l(iat, nn, ncx, ll1, ll2, ll3, l1, l2, l3, myspace, &
2630 rxyz, icell, lstb(iam*myspace + 1), lay, &
2631 rel(1, iam*myspace + 1), cut2, indlst)
2632 lsta(2, iat) = iam*myspace + indlst
2633! write(*,'(a,4(x,i3),100(x,i2))') &
2634! 'iam,iat,lsta',iam,iat,lsta(1,iat),lsta(2,iat), &
2635! (lstb(j),j=lsta(1,iat),lsta(2,iat))
2636 END IF
2637 END DO loop_ii
2638 END DO loop_l1
2639 END DO loop_l2
2640 END DO loop_l3
2641
2642!$OMP ATOMIC UPDATE
2643 indlstx = max(indlstx, indlst)
2644!$OMP END ATOMIC
2645!$OMP END PARALLEL
2646
2647 IF (indlstx >= myspaceout) THEN
2648 WRITE (10, *) count, 'NNBRX too small', nnbrx
2649 DEALLOCATE (lstb, rel)
2650 nnbrx = 3*nnbrx/2
2651 cycle loop_nnbrx
2652 END IF
2653 EXIT loop_nnbrx
2654 END DO loop_nnbrx
2655
2656 istopg = 0
2657!$OMP PARALLEL DEFAULT(NONE) &
2658!$OMP PRIVATE(iam,npr,iat,iat1,iat2,lot,istop,tcoord,tcoord2, &
2659!$OMP tener,tener2,txyz,f2ij,f3ij,f3ik,npjx,npjkx) &
2660!$OMP SHARED (nat,nnbrx,lsta,lstb,rel,ener,ener2,fxyz,coord,coord2,istopg)
2661
2662 npr = 1
2663!$ npr = omp_get_num_threads()
2664 iam = 0
2665!$ iam = omp_get_thread_num()
2666
2667 npjx = 300; npjkx = 6000
2668
2669 IF (npr /= 1) THEN
2670! PARALLEL CASE
2671! create temporary private scalars for reduction sum on energies and
2672! temporary private array for reduction sum on forces
2673!$OMP CRITICAL(omp_eip_lenosky_silicon)
2674 ALLOCATE (txyz(3, nat), f2ij(3, npjx), f3ij(3, npjkx), f3ik(3, npjkx))
2675!$OMP END CRITICAL(omp_eip_lenosky_silicon)
2676 IF (iam == 0) THEN
2677 ener = 0.e0_dp
2678 ener2 = 0.e0_dp
2679 coord = 0.e0_dp
2680 coord2 = 0.e0_dp
2681 END IF
2682!$OMP DO
2683 DO iat = 1, nat
2684 fxyz(1, iat) = 0.e0_dp
2685 fxyz(2, iat) = 0.e0_dp
2686 fxyz(3, iat) = 0.e0_dp
2687 END DO
2688!$OMP BARRIER
2689
2690! Each thread treats at most lot atoms
2691 lot = int(real(nat, kind=dp)/real(npr, kind=dp) + .999999999999e0_dp)
2692 iat1 = iam*lot + 1
2693 iat2 = min((iam + 1)*lot, nat)
2694! write(*,*) 'subfeniat:iat1,iat2,iam',iat1,iat2,iam
2695 CALL subfeniat_l(iat1, iat2, nat, lsta, lstb, rel, tener, tener2, &
2696 tcoord, tcoord2, nnbrx, txyz, f2ij, npjx, f3ij, npjkx, f3ik, istop)
2697!$OMP CRITICAL(omp_eip_lenosky_silicon)
2698 ener = ener + tener
2699 ener2 = ener2 + tener2
2700 coord = coord + tcoord
2701 coord2 = coord2 + tcoord2
2702 istopg = istopg + istop
2703 DO iat = 1, nat
2704 fxyz(1, iat) = fxyz(1, iat) + txyz(1, iat)
2705 fxyz(2, iat) = fxyz(2, iat) + txyz(2, iat)
2706 fxyz(3, iat) = fxyz(3, iat) + txyz(3, iat)
2707 END DO
2708 DEALLOCATE (txyz, f2ij, f3ij, f3ik)
2709!$OMP END CRITICAL(omp_eip_lenosky_silicon)
2710
2711 ELSE
2712! SERIAL CASE
2713 iat1 = 1
2714 iat2 = nat
2715 ALLOCATE (f2ij(3, npjx), f3ij(3, npjkx), f3ik(3, npjkx))
2716 CALL subfeniat_l(iat1, iat2, nat, lsta, lstb, rel, ener, ener2, &
2717 coord, coord2, nnbrx, fxyz, f2ij, npjx, f3ij, npjkx, f3ik, istopg)
2718 DEALLOCATE (f2ij, f3ij, f3ik)
2719
2720 END IF
2721!$OMP END PARALLEL
2722
2723! write(*,*) 'ener,norm force', &
2724! ener,DNRM2(3*nat,fxyz,1)
2725 IF (istopg > 0) cpabort("DIMENSION ERROR (see WARNING above)")
2726 ener_var = ener2/nat - (ener/nat)**2
2727 coord = coord/nat
2728 coord_var = coord2/nat - coord**2
2729
2730 DEALLOCATE (rxyz, icell, lay, lsta, lstb, rel)
2731
2732 END SUBROUTINE eip_lenosky_silicon
2733
2734! **************************************************************************************************
2735!> \brief ...
2736!> \param iat1 ...
2737!> \param iat2 ...
2738!> \param nat ...
2739!> \param lsta ...
2740!> \param lstb ...
2741!> \param rel ...
2742!> \param tener ...
2743!> \param tener2 ...
2744!> \param tcoord ...
2745!> \param tcoord2 ...
2746!> \param nnbrx ...
2747!> \param txyz ...
2748!> \param f2ij ...
2749!> \param npjx ...
2750!> \param f3ij ...
2751!> \param npjkx ...
2752!> \param f3ik ...
2753!> \param istop ...
2754! **************************************************************************************************
2755 SUBROUTINE subfeniat_l(iat1, iat2, nat, lsta, lstb, rel, tener, tener2, &
2756 tcoord, tcoord2, nnbrx, txyz, f2ij, npjx, f3ij, npjkx, f3ik, istop)
2757! for a subset of atoms iat1 to iat2 the routine calculates the (partial) forces
2758! txyz acting on these atoms as well as on the atoms (jat, kat) interacting
2759! with them and their contribution to the energy (tener).
2760! In addition the coordination number tcoord and the second moment of the
2761! local energy tener2 and coordination number tcoord2 are returned
2762 INTEGER :: iat1, iat2, nat, lsta, lstb
2763 REAL(kind=dp) :: rel, tener, tener2, tcoord, tcoord2
2764 INTEGER :: nnbrx
2765 REAL(kind=dp) :: txyz, f2ij
2766 INTEGER :: npjx
2767 REAL(kind=dp) :: f3ij
2768 INTEGER :: npjkx
2769 REAL(kind=dp) :: f3ik
2770 INTEGER :: istop
2771
2772 dimension lsta(2, nat), lstb(nnbrx*nat), rel(5, nnbrx*nat), txyz(3, nat)
2773 dimension f2ij(3, npjx), f3ij(3, npjkx), f3ik(3, npjkx)
2774 REAL(kind=dp), PARAMETER :: tmin_phi = 0.1500000e+01_dp
2775 REAL(kind=dp), PARAMETER :: tmax_phi = 0.4500000e+01_dp
2776 REAL(kind=dp), PARAMETER :: hi_phi = 3.00000000000e0_dp
2777 REAL(kind=dp), PARAMETER :: hsixth_phi = 5.55555555555556e-002_dp
2778 REAL(kind=dp), PARAMETER :: h2sixth_phi = 1.85185185185185e-002_dp
2779 REAL(kind=dp), PARAMETER, DIMENSION(0:9) :: cof_phi = &
2780 [0.69299400000000e+01_dp, -0.43995000000000e+00_dp, &
2781 -0.17012300000000e+01_dp, -0.16247300000000e+01_dp, &
2782 -0.99696000000000e+00_dp, -0.27391000000000e+00_dp, &
2783 -0.24990000000000e-01_dp, -0.17840000000000e-01_dp, &
2784 -0.96100000000000e-02_dp, 0.00000000000000e+00_dp]
2785 REAL(kind=dp), PARAMETER, DIMENSION(0:9) :: dof_phi = &
2786 [0.16533229480429e+03_dp, 0.39415410391417e+02_dp, &
2787 0.68710036300407e+01_dp, 0.53406950884203e+01_dp, &
2788 0.15347960162782e+01_dp, -0.63347591535331e+01_dp, &
2789 -0.17987794021458e+01_dp, 0.47429676211617e+00_dp, &
2790 -0.40087646318907e-01_dp, -0.23942617684055e+00_dp]
2791 REAL(kind=dp), PARAMETER :: tmin_rho = 0.1500000e+01_dp
2792 REAL(kind=dp), PARAMETER :: tmax_rho = 0.3500000e+01_dp
2793 REAL(kind=dp), PARAMETER :: hi_rho = 5.00000000000e0_dp
2794 REAL(kind=dp), PARAMETER :: hsixth_rho = 3.33333333333333e-002_dp
2795 REAL(kind=dp), PARAMETER :: h2sixth_rho = 6.66666666666667e-003_dp
2796 REAL(kind=dp), PARAMETER, DIMENSION(0:10) :: cof_rho = &
2797 [0.13747000000000e+00_dp, -0.14831000000000e+00_dp, &
2798 -0.55972000000000e+00_dp, -0.73110000000000e+00_dp, &
2799 -0.76283000000000e+00_dp, -0.72918000000000e+00_dp, &
2800 -0.66620000000000e+00_dp, -0.57328000000000e+00_dp, &
2801 -0.40690000000000e+00_dp, -0.16662000000000e+00_dp, &
2802 0.00000000000000e+00_dp]
2803 REAL(kind=dp), PARAMETER, DIMENSION(0:10) :: dof_rho = &
2804 [-0.32275496741918e+01_dp, -0.64119006516165e+01_dp, &
2805 0.10030652280658e+02_dp, 0.22937915289857e+01_dp, &
2806 0.17416816033995e+01_dp, 0.54648205741626e+00_dp, &
2807 0.47189016693543e+00_dp, 0.20569572748420e+01_dp, &
2808 0.23192807336964e+01_dp, -0.24908020962757e+00_dp, &
2809 -0.12371959895186e+02_dp]
2810 REAL(kind=dp), PARAMETER :: tmin_fff = 0.1500000e+01_dp
2811 REAL(kind=dp), PARAMETER :: tmax_fff = 0.3500000e+01_dp
2812 REAL(kind=dp), PARAMETER :: hi_fff = 4.50000000000e0_dp
2813 REAL(kind=dp), PARAMETER :: hsixth_fff = 3.70370370370370e-002_dp
2814 REAL(kind=dp), PARAMETER :: h2sixth_fff = 8.23045267489712e-003_dp
2815 REAL(kind=dp), PARAMETER, DIMENSION(0:9) :: cof_fff = &
2816 [0.12503100000000e+01_dp, 0.86821000000000e+00_dp, &
2817 0.60846000000000e+00_dp, 0.48756000000000e+00_dp, &
2818 0.44163000000000e+00_dp, 0.37610000000000e+00_dp, &
2819 0.27145000000000e+00_dp, 0.14814000000000e+00_dp, &
2820 0.48550000000000e-01_dp, 0.00000000000000e+00_dp]
2821 REAL(kind=dp), PARAMETER, DIMENSION(0:9) :: dof_fff = &
2822 [0.27904652711432e+02_dp, -0.45230754228635e+01_dp, &
2823 0.50531739800222e+01_dp, 0.11806545027747e+01_dp, &
2824 -0.66693699112098e+00_dp, -0.89430653829079e+00_dp, &
2825 -0.50891685571587e+00_dp, 0.66278396115427e+00_dp, &
2826 0.73976101109878e+00_dp, 0.25795319944506e+01_dp]
2827 REAL(kind=dp), PARAMETER :: tmin_uuu = -0.1770930e+01_dp
2828 REAL(kind=dp), PARAMETER :: tmax_uuu = 0.7908520e+01_dp
2829 REAL(kind=dp), PARAMETER :: hi_uuu = 0.723181585730594e0_dp
2830 REAL(kind=dp), PARAMETER :: hsixth_uuu = 0.230463095238095e0_dp
2831 REAL(kind=dp), PARAMETER :: h2sixth_uuu = 0.318679429600340e0_dp
2832 REAL(kind=dp), PARAMETER, DIMENSION(0:7) :: cof_uuu = &
2833 [-0.10749300000000e+01_dp, -0.20045000000000e+00_dp, &
2834 0.41422000000000e+00_dp, 0.87939000000000e+00_dp, &
2835 0.12668900000000e+01_dp, 0.16299800000000e+01_dp, &
2836 0.19773800000000e+01_dp, 0.23961800000000e+01_dp]
2837 REAL(kind=dp), PARAMETER, DIMENSION(0:7) :: dof_uuu = &
2838 [-0.14827125747284e+00_dp, -0.14922155328475e+00_dp, &
2839 -0.70113224223509e-01_dp, -0.39449020349230e-01_dp, &
2840 -0.15815242579643e-01_dp, 0.26112640061855e-01_dp, &
2841 -0.13786974745095e+00_dp, 0.74941595372657e+00_dp]
2842 REAL(kind=dp), PARAMETER :: tmin_ggg = -0.1000000e+01_dp
2843 REAL(kind=dp), PARAMETER :: tmax_ggg = 0.8001400e+00_dp
2844 REAL(kind=dp), PARAMETER :: hi_ggg = 3.88858644327663e0_dp
2845 REAL(kind=dp), PARAMETER :: hsixth_ggg = 4.28604761904762e-002_dp
2846 REAL(kind=dp), PARAMETER :: h2sixth_ggg = 1.10221225156463e-002_dp
2847 REAL(kind=dp), PARAMETER, DIMENSION(0:7) :: cof_ggg = &
2848 [0.52541600000000e+01_dp, 0.23591500000000e+01_dp, &
2849 0.11959500000000e+01_dp, 0.12299500000000e+01_dp, &
2850 0.20356500000000e+01_dp, 0.34247400000000e+01_dp, &
2851 0.49485900000000e+01_dp, 0.56179900000000e+01_dp]
2852 REAL(kind=dp), PARAMETER, DIMENSION(0:7) :: dof_ggg = &
2853 [0.15826876132396e+02_dp, 0.31176239377907e+02_dp, &
2854 0.16589446539683e+02_dp, 0.11083892500520e+02_dp, &
2855 0.90887216383860e+01_dp, 0.54902279653967e+01_dp, &
2856 -0.18823313223755e+02_dp, -0.77183416481005e+01_dp]
2857
2858 REAL(kind=dp) :: a2_fff, a2_ggg, a_fff, a_ggg, b2_fff, b2_ggg, b_fff, &
2859 b_ggg, cof1_fff, cof1_ggg, cof2_fff, cof2_ggg, cof3_fff, &
2860 cof3_ggg, cof4_fff, cof4_ggg, cof_fff_khi, cof_fff_klo, &
2861 cof_ggg_khi, cof_ggg_klo, coord_iat, costheta, dens, &
2862 dens2, dens3, dof_fff_khi, dof_fff_klo, dof_ggg_khi, &
2863 dof_ggg_klo, e_phi, e_uuu, ener_iat, ep_phi, ep_uuu, &
2864 fij, fijp, fik, fikp, fxij, fxik, fyij, fyik, fzij, fzik, &
2865 gjik, gjikp, rho, rhop, rij, rik, sij, sik, t1, t2, t3, t4, &
2866 tt, tt_fff, tt_ggg, xarg, ypt1_fff, ypt1_ggg, ypt2_fff, &
2867 ypt2_ggg, yt1_fff, yt1_ggg, yt2_fff, yt2_ggg
2868
2869 INTEGER :: iat, jat, jbr, jcnt, jkcnt, kat, kbr, khi_fff, khi_ggg, &
2870 klo_fff, klo_ggg
2871
2872! initialize temporary private scalars for reduction sum on energies and
2873! private workarray txyz for forces forces
2874 tener = 0.e0_dp
2875 tener2 = 0.e0_dp
2876 tcoord = 0.e0_dp
2877 tcoord2 = 0.e0_dp
2878 istop = 0
2879 DO iat = 1, nat
2880 txyz(1, iat) = 0.e0_dp
2881 txyz(2, iat) = 0.e0_dp
2882 txyz(3, iat) = 0.e0_dp
2883 END DO
2884
2885! calculation of forces, energy
2886
2887 forces_and_energy: DO iat = iat1, iat2
2888
2889 dens2 = 0.e0_dp
2890 dens3 = 0.e0_dp
2891 jcnt = 0
2892 jkcnt = 0
2893 coord_iat = 0.e0_dp
2894 ener_iat = 0.e0_dp
2895 calculate: DO jbr = lsta(1, iat), lsta(2, iat)
2896 jat = lstb(jbr)
2897 jcnt = jcnt + 1
2898 IF (jcnt > npjx) THEN
2899 WRITE (*, *) 'WARNING: enlarge npjx'
2900 istop = 1
2901 RETURN
2902 END IF
2903
2904 fxij = rel(1, jbr)
2905 fyij = rel(2, jbr)
2906 fzij = rel(3, jbr)
2907 rij = rel(4, jbr)
2908 sij = rel(5, jbr)
2909
2910! coordination number calculated with soft cutoff between first
2911! nearest neighbor and midpoint of first and second nearest neighbor
2912 IF (rij <= 2.36e0_dp) THEN
2913 coord_iat = coord_iat + 1.e0_dp
2914 ELSE IF (rij >= 3.12e0_dp) THEN
2915 ELSE
2916 xarg = (rij - 2.36e0_dp)*(1.e0_dp/(3.12e0_dp - 2.36e0_dp))
2917 coord_iat = coord_iat + (2*xarg + 1.e0_dp)*(xarg - 1.e0_dp)**2
2918 END IF
2919
2920! pairpotential term
2921 CALL splint(cof_phi, dof_phi, tmin_phi, tmax_phi, &
2922 hsixth_phi, h2sixth_phi, hi_phi, 10, rij, e_phi, ep_phi)
2923 ener_iat = ener_iat + (e_phi*.5e0_dp)
2924 txyz(1, iat) = txyz(1, iat) - fxij*(ep_phi*.5e0_dp)
2925 txyz(2, iat) = txyz(2, iat) - fyij*(ep_phi*.5e0_dp)
2926 txyz(3, iat) = txyz(3, iat) - fzij*(ep_phi*.5e0_dp)
2927 txyz(1, jat) = txyz(1, jat) + fxij*(ep_phi*.5e0_dp)
2928 txyz(2, jat) = txyz(2, jat) + fyij*(ep_phi*.5e0_dp)
2929 txyz(3, jat) = txyz(3, jat) + fzij*(ep_phi*.5e0_dp)
2930
2931! 2 body embedding term
2932 CALL splint(cof_rho, dof_rho, tmin_rho, tmax_rho, &
2933 hsixth_rho, h2sixth_rho, hi_rho, 11, rij, rho, rhop)
2934 dens2 = dens2 + rho
2935 f2ij(1, jcnt) = fxij*rhop
2936 f2ij(2, jcnt) = fyij*rhop
2937 f2ij(3, jcnt) = fzij*rhop
2938
2939! 3 body embedding term
2940 CALL splint(cof_fff, dof_fff, tmin_fff, tmax_fff, &
2941 hsixth_fff, h2sixth_fff, hi_fff, 10, rij, fij, fijp)
2942
2943 embed_3body: DO kbr = lsta(1, iat), lsta(2, iat)
2944 kat = lstb(kbr)
2945 IF (kat < jat) THEN
2946 jkcnt = jkcnt + 1
2947 IF (jkcnt > npjkx) THEN
2948 WRITE (*, *) 'WARNING: enlarge npjkx', npjkx
2949 istop = 1
2950 RETURN
2951 END IF
2952
2953! begin unoptimized original version:
2954! fxik=rel(1,kbr)
2955! fyik=rel(2,kbr)
2956! fzik=rel(3,kbr)
2957! rik=rel(4,kbr)
2958! sik=rel(5,kbr)
2959!
2960! call splint(cof_fff,dof_fff,tmin_fff,tmax_fff, &
2961! hsixth_fff,h2sixth_fff,hi_fff,10,rik,fik,fikp)
2962! costheta=fxij*fxik+fyij*fyik+fzij*fzik
2963! call splint(cof_ggg,dof_ggg,tmin_ggg,tmax_ggg, &
2964! hsixth_ggg,h2sixth_ggg,hi_ggg,8,costheta,gjik,gjikp)
2965! end unoptimized original version:
2966
2967! begin optimized version
2968 rik = rel(4, kbr)
2969 IF (rik > tmax_fff) THEN
2970 fikp = 0.e0_dp; fik = 0.e0_dp
2971 gjik = 0.e0_dp; gjikp = 0.e0_dp; sik = 0.e0_dp
2972 costheta = 0.e0_dp; fxik = 0.e0_dp; fyik = 0.e0_dp; fzik = 0.e0_dp
2973 ELSE IF (rik < tmin_fff) THEN
2974 fxik = rel(1, kbr)
2975 fyik = rel(2, kbr)
2976 fzik = rel(3, kbr)
2977 costheta = fxij*fxik + fyij*fyik + fzij*fzik
2978 sik = rel(5, kbr)
2979 fikp = hi_fff*(cof_fff(1) - cof_fff(0)) - &
2980 (dof_fff(1) + 2.e0_dp*dof_fff(0))*hsixth_fff
2981 fik = cof_fff(0) + (rik - tmin_fff)*fikp
2982 tt_ggg = (costheta - tmin_ggg)*hi_ggg
2983 IF (costheta > tmax_ggg) THEN
2984 gjikp = hi_ggg*(cof_ggg(8 - 1) - cof_ggg(8 - 2)) + &
2985 (2.e0_dp*dof_ggg(8 - 1) + dof_ggg(8 - 2))*hsixth_ggg
2986 gjik = cof_ggg(8 - 1) + (costheta - tmax_ggg)*gjikp
2987 ELSE
2988 klo_ggg = int(tt_ggg)
2989 khi_ggg = klo_ggg + 1
2990 cof_ggg_klo = cof_ggg(klo_ggg)
2991 dof_ggg_klo = dof_ggg(klo_ggg)
2992 b_ggg = tt_ggg - klo_ggg
2993 a_ggg = 1.e0_dp - b_ggg
2994 cof_ggg_khi = cof_ggg(khi_ggg)
2995 dof_ggg_khi = dof_ggg(khi_ggg)
2996 b2_ggg = b_ggg*b_ggg
2997 gjik = a_ggg*cof_ggg_klo
2998 gjikp = cof_ggg_khi - cof_ggg_klo
2999 a2_ggg = a_ggg*a_ggg
3000 cof1_ggg = a2_ggg - 1.e0_dp
3001 cof2_ggg = b2_ggg - 1.e0_dp
3002 gjik = gjik + b_ggg*cof_ggg_khi
3003 gjikp = hi_ggg*gjikp
3004 cof3_ggg = 3.e0_dp*b2_ggg
3005 cof4_ggg = 3.e0_dp*a2_ggg
3006 cof1_ggg = a_ggg*cof1_ggg
3007 cof2_ggg = b_ggg*cof2_ggg
3008 cof3_ggg = cof3_ggg - 1.e0_dp
3009 cof4_ggg = cof4_ggg - 1.e0_dp
3010 yt1_ggg = cof1_ggg*dof_ggg_klo
3011 yt2_ggg = cof2_ggg*dof_ggg_khi
3012 ypt1_ggg = cof3_ggg*dof_ggg_khi
3013 ypt2_ggg = cof4_ggg*dof_ggg_klo
3014 gjik = gjik + (yt1_ggg + yt2_ggg)*h2sixth_ggg
3015 gjikp = gjikp + (ypt1_ggg - ypt2_ggg)*hsixth_ggg
3016 END IF
3017 ELSE
3018 fxik = rel(1, kbr)
3019 tt_fff = rik - tmin_fff
3020 costheta = fxij*fxik
3021 fyik = rel(2, kbr)
3022 tt_fff = tt_fff*hi_fff
3023 costheta = costheta + fyij*fyik
3024 fzik = rel(3, kbr)
3025 klo_fff = int(tt_fff)
3026 costheta = costheta + fzij*fzik
3027 sik = rel(5, kbr)
3028 tt_ggg = (costheta - tmin_ggg)*hi_ggg
3029 IF (costheta > tmax_ggg) THEN
3030 gjikp = hi_ggg*(cof_ggg(8 - 1) - cof_ggg(8 - 2)) + &
3031 (2.e0_dp*dof_ggg(8 - 1) + dof_ggg(8 - 2))*hsixth_ggg
3032 gjik = cof_ggg(8 - 1) + (costheta - tmax_ggg)*gjikp
3033 khi_fff = klo_fff + 1
3034 cof_fff_klo = cof_fff(klo_fff)
3035 dof_fff_klo = dof_fff(klo_fff)
3036 b_fff = tt_fff - klo_fff
3037 a_fff = 1.e0_dp - b_fff
3038 cof_fff_khi = cof_fff(khi_fff)
3039 dof_fff_khi = dof_fff(khi_fff)
3040 b2_fff = b_fff*b_fff
3041 fik = a_fff*cof_fff_klo
3042 fikp = cof_fff_khi - cof_fff_klo
3043 a2_fff = a_fff*a_fff
3044 cof1_fff = a2_fff - 1.e0_dp
3045 cof2_fff = b2_fff - 1.e0_dp
3046 fik = fik + b_fff*cof_fff_khi
3047 fikp = hi_fff*fikp
3048 cof3_fff = 3.e0_dp*b2_fff
3049 cof4_fff = 3.e0_dp*a2_fff
3050 cof1_fff = a_fff*cof1_fff
3051 cof2_fff = b_fff*cof2_fff
3052 cof3_fff = cof3_fff - 1.e0_dp
3053 cof4_fff = cof4_fff - 1.e0_dp
3054 yt1_fff = cof1_fff*dof_fff_klo
3055 yt2_fff = cof2_fff*dof_fff_khi
3056 ypt1_fff = cof3_fff*dof_fff_khi
3057 ypt2_fff = cof4_fff*dof_fff_klo
3058 fik = fik + (yt1_fff + yt2_fff)*h2sixth_fff
3059 fikp = fikp + (ypt1_fff - ypt2_fff)*hsixth_fff
3060 ELSE
3061 klo_ggg = int(tt_ggg)
3062 khi_ggg = klo_ggg + 1
3063 khi_fff = klo_fff + 1
3064 cof_ggg_klo = cof_ggg(klo_ggg)
3065 cof_fff_klo = cof_fff(klo_fff)
3066 dof_ggg_klo = dof_ggg(klo_ggg)
3067 dof_fff_klo = dof_fff(klo_fff)
3068 b_ggg = tt_ggg - klo_ggg
3069 b_fff = tt_fff - klo_fff
3070 a_ggg = 1.e0_dp - b_ggg
3071 a_fff = 1.e0_dp - b_fff
3072 cof_ggg_khi = cof_ggg(khi_ggg)
3073 cof_fff_khi = cof_fff(khi_fff)
3074 dof_ggg_khi = dof_ggg(khi_ggg)
3075 dof_fff_khi = dof_fff(khi_fff)
3076 b2_ggg = b_ggg*b_ggg
3077 b2_fff = b_fff*b_fff
3078 gjik = a_ggg*cof_ggg_klo
3079 fik = a_fff*cof_fff_klo
3080 gjikp = cof_ggg_khi - cof_ggg_klo
3081 fikp = cof_fff_khi - cof_fff_klo
3082 a2_ggg = a_ggg*a_ggg
3083 a2_fff = a_fff*a_fff
3084 cof1_ggg = a2_ggg - 1.e0_dp
3085 cof1_fff = a2_fff - 1.e0_dp
3086 cof2_ggg = b2_ggg - 1.e0_dp
3087 cof2_fff = b2_fff - 1.e0_dp
3088 gjik = gjik + b_ggg*cof_ggg_khi
3089 fik = fik + b_fff*cof_fff_khi
3090 gjikp = hi_ggg*gjikp
3091 fikp = hi_fff*fikp
3092 cof3_ggg = 3.e0_dp*b2_ggg
3093 cof3_fff = 3.e0_dp*b2_fff
3094 cof4_ggg = 3.e0_dp*a2_ggg
3095 cof4_fff = 3.e0_dp*a2_fff
3096 cof1_ggg = a_ggg*cof1_ggg
3097 cof1_fff = a_fff*cof1_fff
3098 cof2_ggg = b_ggg*cof2_ggg
3099 cof2_fff = b_fff*cof2_fff
3100 cof3_ggg = cof3_ggg - 1.e0_dp
3101 cof3_fff = cof3_fff - 1.e0_dp
3102 cof4_ggg = cof4_ggg - 1.e0_dp
3103 cof4_fff = cof4_fff - 1.e0_dp
3104 yt1_ggg = cof1_ggg*dof_ggg_klo
3105 yt1_fff = cof1_fff*dof_fff_klo
3106 yt2_ggg = cof2_ggg*dof_ggg_khi
3107 yt2_fff = cof2_fff*dof_fff_khi
3108 ypt1_ggg = cof3_ggg*dof_ggg_khi
3109 ypt1_fff = cof3_fff*dof_fff_khi
3110 ypt2_ggg = cof4_ggg*dof_ggg_klo
3111 ypt2_fff = cof4_fff*dof_fff_klo
3112 gjik = gjik + (yt1_ggg + yt2_ggg)*h2sixth_ggg
3113 fik = fik + (yt1_fff + yt2_fff)*h2sixth_fff
3114 gjikp = gjikp + (ypt1_ggg - ypt2_ggg)*hsixth_ggg
3115 fikp = fikp + (ypt1_fff - ypt2_fff)*hsixth_fff
3116 END IF
3117 END IF
3118! end optimized version
3119
3120 tt = fij*fik
3121 dens3 = dens3 + tt*gjik
3122
3123 t1 = fijp*fik*gjik
3124 t2 = sij*(tt*gjikp)
3125 f3ij(1, jkcnt) = fxij*t1 + (fxik - fxij*costheta)*t2
3126 f3ij(2, jkcnt) = fyij*t1 + (fyik - fyij*costheta)*t2
3127 f3ij(3, jkcnt) = fzij*t1 + (fzik - fzij*costheta)*t2
3128
3129 t3 = fikp*fij*gjik
3130 t4 = sik*(tt*gjikp)
3131 f3ik(1, jkcnt) = fxik*t3 + (fxij - fxik*costheta)*t4
3132 f3ik(2, jkcnt) = fyik*t3 + (fyij - fyik*costheta)*t4
3133 f3ik(3, jkcnt) = fzik*t3 + (fzij - fzik*costheta)*t4
3134 END IF
3135
3136 END DO embed_3body
3137 END DO calculate
3138
3139 dens = dens2 + dens3
3140 CALL splint(cof_uuu, dof_uuu, tmin_uuu, tmax_uuu, &
3141 hsixth_uuu, h2sixth_uuu, hi_uuu, 8, dens, e_uuu, ep_uuu)
3142 ener_iat = ener_iat + e_uuu
3143
3144! Only now ep_uu is known and the forces can be calculated, lets loop again
3145 jcnt = 0
3146 jkcnt = 0
3147 loop_again: DO jbr = lsta(1, iat), lsta(2, iat)
3148 jat = lstb(jbr)
3149 jcnt = jcnt + 1
3150 txyz(1, iat) = txyz(1, iat) - ep_uuu*f2ij(1, jcnt)
3151 txyz(2, iat) = txyz(2, iat) - ep_uuu*f2ij(2, jcnt)
3152 txyz(3, iat) = txyz(3, iat) - ep_uuu*f2ij(3, jcnt)
3153 txyz(1, jat) = txyz(1, jat) + ep_uuu*f2ij(1, jcnt)
3154 txyz(2, jat) = txyz(2, jat) + ep_uuu*f2ij(2, jcnt)
3155 txyz(3, jat) = txyz(3, jat) + ep_uuu*f2ij(3, jcnt)
3156
3157! 3 body embedding term
3158 DO kbr = lsta(1, iat), lsta(2, iat)
3159 kat = lstb(kbr)
3160 IF (kat < jat) THEN
3161 jkcnt = jkcnt + 1
3162
3163 txyz(1, iat) = txyz(1, iat) - ep_uuu*(f3ij(1, jkcnt) + f3ik(1, jkcnt))
3164 txyz(2, iat) = txyz(2, iat) - ep_uuu*(f3ij(2, jkcnt) + f3ik(2, jkcnt))
3165 txyz(3, iat) = txyz(3, iat) - ep_uuu*(f3ij(3, jkcnt) + f3ik(3, jkcnt))
3166 txyz(1, jat) = txyz(1, jat) + ep_uuu*f3ij(1, jkcnt)
3167 txyz(2, jat) = txyz(2, jat) + ep_uuu*f3ij(2, jkcnt)
3168 txyz(3, jat) = txyz(3, jat) + ep_uuu*f3ij(3, jkcnt)
3169 txyz(1, kat) = txyz(1, kat) + ep_uuu*f3ik(1, jkcnt)
3170 txyz(2, kat) = txyz(2, kat) + ep_uuu*f3ik(2, jkcnt)
3171 txyz(3, kat) = txyz(3, kat) + ep_uuu*f3ik(3, jkcnt)
3172 END IF
3173 END DO
3174
3175 END DO loop_again
3176
3177! write(*,'(a,i4,x,e19.12,x,e10.3)') 'iat,ener_iat,coord_iat', &
3178! iat,ener_iat,coord_iat
3179 tener = tener + ener_iat
3180 tener2 = tener2 + ener_iat**2
3181 tcoord = tcoord + coord_iat
3182 tcoord2 = tcoord2 + coord_iat**2
3183
3184 END DO forces_and_energy
3185
3186 END SUBROUTINE subfeniat_l
3187
3188! **************************************************************************************************
3189!> \brief ...
3190!> \param iat ...
3191!> \param nn ...
3192!> \param ncx ...
3193!> \param ll1 ...
3194!> \param ll2 ...
3195!> \param ll3 ...
3196!> \param l1 ...
3197!> \param l2 ...
3198!> \param l3 ...
3199!> \param myspace ...
3200!> \param rxyz ...
3201!> \param icell ...
3202!> \param lstb ...
3203!> \param lay ...
3204!> \param rel ...
3205!> \param cut2 ...
3206!> \param indlst ...
3207! **************************************************************************************************
3208 SUBROUTINE sublstiat_l(iat, nn, ncx, ll1, ll2, ll3, l1, l2, l3, myspace, &
3209 rxyz, icell, lstb, lay, rel, cut2, indlst)
3210! finds the neighbours of atom iat (specified by lsta and lstb) and and
3211! the relative position rel of iat with respect to these neighbours
3212 INTEGER :: iat, nn, ncx, ll1, ll2, ll3, l1, l2, l3, &
3213 myspace
3214 REAL(kind=dp) :: rxyz
3215 INTEGER :: icell, lstb, lay
3216 REAL(kind=dp) :: rel, cut2
3217 INTEGER :: indlst
3218
3219 dimension rxyz(3, nn), lay(nn), icell(0:ncx, -1:ll1, -1:ll2, -1:ll3), &
3220 lstb(0:myspace - 1), rel(5, 0:myspace - 1)
3221
3222 INTEGER :: jat, jj, k1, k2, k3
3223 REAL(kind=dp) :: rr2, tt, xrel, yrel, zrel, tti
3224
3225 loop_k3: DO k3 = l3 - 1, l3 + 1
3226 loop_k2: DO k2 = l2 - 1, l2 + 1
3227 loop_k1: DO k1 = l1 - 1, l1 + 1
3228 loop_jj: DO jj = 1, icell(0, k1, k2, k3)
3229 jat = icell(jj, k1, k2, k3)
3230 IF (jat == iat) cycle loop_k3
3231 xrel = rxyz(1, iat) - rxyz(1, jat)
3232 yrel = rxyz(2, iat) - rxyz(2, jat)
3233 zrel = rxyz(3, iat) - rxyz(3, jat)
3234 rr2 = xrel**2 + yrel**2 + zrel**2
3235 IF (rr2 <= cut2) THEN
3236 indlst = min(indlst, myspace - 1)
3237 lstb(indlst) = lay(jat)
3238! write(*,*) 'iat,indlst,lay(jat)',iat,indlst,lay(jat)
3239 tt = sqrt(rr2)
3240 tti = 1.e0_dp/tt
3241 rel(1, indlst) = xrel*tti
3242 rel(2, indlst) = yrel*tti
3243 rel(3, indlst) = zrel*tti
3244 rel(4, indlst) = tt
3245 rel(5, indlst) = tti
3246 indlst = indlst + 1
3247 END IF
3248 END DO loop_jj
3249 END DO loop_k1
3250 END DO loop_k2
3251 END DO loop_k3
3252
3253 RETURN
3254 END SUBROUTINE sublstiat_l
3255
3256! **************************************************************************************************
3257!> \brief ...
3258!> \param ya ...
3259!> \param y2a ...
3260!> \param tmin ...
3261!> \param tmax ...
3262!> \param hsixth ...
3263!> \param h2sixth ...
3264!> \param hi ...
3265!> \param n ...
3266!> \param x ...
3267!> \param y ...
3268!> \param yp ...
3269! **************************************************************************************************
3270 SUBROUTINE splint(ya, y2a, tmin, tmax, hsixth, h2sixth, hi, n, x, y, yp)
3271 REAL(kind=dp) :: ya, y2a, tmin, tmax, hsixth, h2sixth, hi
3272 INTEGER :: n
3273 REAL(kind=dp) :: x, y, yp
3274
3275 dimension y2a(0:n - 1), ya(0:n - 1)
3276 REAL(kind=dp) :: a, a2, b, b2, cof1, cof2, cof3, cof4, tt, &
3277 y2a_khi, ya_klo, y2a_klo, ya_khi, ypt1, ypt2, yt1, yt2
3278 INTEGER :: klo, khi
3279
3280! interpolate if the argument is outside the cubic spline interval [tmin,tmax]
3281 tt = (x - tmin)*hi
3282 IF (x < tmin) THEN
3283 yp = hi*(ya(1) - ya(0)) - &
3284 (y2a(1) + 2.e0_dp*y2a(0))*hsixth
3285 y = ya(0) + (x - tmin)*yp
3286 ELSE IF (x > tmax) THEN
3287 yp = hi*(ya(n - 1) - ya(n - 2)) + &
3288 (2.e0_dp*y2a(n - 1) + y2a(n - 2))*hsixth
3289 y = ya(n - 1) + (x - tmax)*yp
3290! otherwise evaluate cubic spline
3291 ELSE
3292 klo = int(tt)
3293 khi = klo + 1
3294 ya_klo = ya(klo)
3295 y2a_klo = y2a(klo)
3296 b = tt - klo
3297 a = 1.e0_dp - b
3298 ya_khi = ya(khi)
3299 y2a_khi = y2a(khi)
3300 b2 = b*b
3301 y = a*ya_klo
3302 yp = ya_khi - ya_klo
3303 a2 = a*a
3304 cof1 = a2 - 1.e0_dp
3305 cof2 = b2 - 1.e0_dp
3306 y = y + b*ya_khi
3307 yp = hi*yp
3308 cof3 = 3.e0_dp*b2
3309 cof4 = 3.e0_dp*a2
3310 cof1 = a*cof1
3311 cof2 = b*cof2
3312 cof3 = cof3 - 1.e0_dp
3313 cof4 = cof4 - 1.e0_dp
3314 yt1 = cof1*y2a_klo
3315 yt2 = cof2*y2a_khi
3316 ypt1 = cof3*y2a_khi
3317 ypt2 = cof4*y2a_klo
3318 y = y + (yt1 + yt2)*h2sixth
3319 yp = yp + (ypt1 - ypt2)*hsixth
3320 END IF
3321 RETURN
3322 END SUBROUTINE splint
3323
3324! **************************************************************************************************
3325! Additional EIP kernels consolidated here to match the original Bazant/Lenosky layout.
3326! **************************************************************************************************
3327
3328! **************************************************************************************************
3329!> \brief ...
3330!> \param nat ...
3331!> \param alat ...
3332!> \param rxyz0 ...
3333!> \param fxyz ...
3334!> \param etot ...
3335!> \param count ...
3336! **************************************************************************************************
3337 SUBROUTINE eip_stillinger_weber_silicon(nat, alat, rxyz0, fxyz, etot, count)
3338!*****************************************************************************************
3339! This subroutine evaluates the Stillinger Weber Silicon potential with linear scaling
3340! COPYRIGHT
3341! Copyright (C) 2009 AIST, UNIBAS
3342! This file is distributed under the terms of the
3343! GNU General Public License, see
3344! http://www.gnu.org/copyleft/gpl.txt .
3345!
3346! Implementation: Original version was written by Tetsuya Morishita, AIST Tsukuba (JP)
3347! Improved by M. Amsler, S. Goedecker, Basel University (CH), 2009
3348!
3349! Note:
3350!
3351! aa is the parameter A given on page 5263 of PRB 31, 5262 (1985).
3352! bb is the parameter B given on page 5263 of PRB 31, 5262 (1985).
3353! ra is the parameter a given on page 5263 of PRB 31, 5262 (1985).
3354! gam and ramda are the parameters gamma and lambda
3355! for the 3-body term, respectively (see Eq. (2.5) in the paper).
3356!
3357! Input:
3358! nat, integer: the number of atoms
3359! alat, real(8), dim(3) : the three edges of the orthoromic simulation cell, periodic boundaries are applied
3360! and atoms outside the cell will be brought back into the box
3361! rxyz, real(8), dim(3,nat) : the xyz cartesian components of the atomic positions in Angstroem
3362!
3363! Output:
3364! fxyz, real(8), dim(3,nat): the xyz cartesian forces in on the corresponding atomic components n eV/A
3365! etot, real(8) : total potential energy, 2-body and 3-body, in eV
3366! count, real(8) : increased by 1.d0 at each call of this subroutine,
3367! needs to be initialized to 0.d0 before calling this routine for the first time
3368!
3369! Other variables:
3370! p: the 2-body potential energy
3371! p3: the 3-body potential energy
3372! fx(i) , fy(i) , fz(i) are the 2-body forces on atom i.
3373! fx3(i), fy3(i), fz3(i) are the 3-body forces on atom i.
3374! fxyz(3,nat) contains both 2-body and 3-body forces
3375!
3376! All units follow the description in PRB 31, 5262 (1985).
3377!*****************************************************************************************
3378
3379 INTEGER :: nat
3380 REAL(8) :: alat(3), rxyz0(3, nat), fxyz(3, nat), &
3381 etot, count
3382
3383 INTEGER :: i, iam, iat, ii, il, in, indlst, &
3384 indlstx, ipb, l1, l2, l3, laymx, ll1, &
3385 ll2, ll3, myspace, myspaceout, ncx, &
3386 ndat, nn, nnbrx, npjkx, npjx, npr
3387 INTEGER, ALLOCATABLE, DIMENSION(:) :: lay, lstb
3388 INTEGER, ALLOCATABLE, DIMENSION(:, :) :: lsta
3389 INTEGER, ALLOCATABLE, DIMENSION(:, :, :, :) :: icell
3390 REAL(8) :: cut, cut2, eps, esigma, fx(nat), &
3391 fx3(nat), fy(nat), fy3(nat), fz(nat), &
3392 fz3(nat), isigma, p, p3, pv3, ra, &
3393 rlc1i, rlc2i, rlc3i, sigma
3394 REAL(8), ALLOCATABLE, DIMENSION(:, :) :: rel, rxyz
3395
3396 parameter(ra=1.8d0)
3397 parameter(sigma=2.0951d0, eps=2.167239428587d0)
3398
3399 count = count + 1.d0
3400 cut = sigma*ra*2.d0
3401 isigma = 1.d0/sigma
3402 esigma = eps*isigma
3403
3404! linear scaling calculation of verlet list, only serial
3405 ll1 = int(alat(1)/cut)
3406 IF (ll1 < 1) cpabort("alat(1) too small")
3407 ll2 = int(alat(2)/cut)
3408 IF (ll2 < 1) cpabort("alat(2) too small")
3409 ll3 = int(alat(3)/cut)
3410 IF (ll3 < 1) cpabort("alat(3) too small")
3411
3412 npr = 1
3413 ncx = 29
3414 DO
3415 ncx = ncx*2
3416 ALLOCATE (icell(0:ncx, -1:ll1, -1:ll2, -1:ll3))
3417 icell(0, :, :, :) = 0
3418 rlc1i = ll1/alat(1)
3419 rlc2i = ll2/alat(2)
3420 rlc3i = ll3/alat(3)
3421
3422 DO iat = 1, nat
3423 rxyz0(1, iat) = modulo(modulo(rxyz0(1, iat), alat(1)), alat(1))
3424 rxyz0(2, iat) = modulo(modulo(rxyz0(2, iat), alat(2)), alat(2))
3425 rxyz0(3, iat) = modulo(modulo(rxyz0(3, iat), alat(3)), alat(3))
3426 l1 = int(rxyz0(1, iat)*rlc1i)
3427 l2 = int(rxyz0(2, iat)*rlc2i)
3428 l3 = int(rxyz0(3, iat)*rlc3i)
3429
3430 ii = icell(0, l1, l2, l3)
3431 ii = ii + 1
3432 icell(0, l1, l2, l3) = ii
3433 IF (ii > ncx) THEN
3434 DEALLOCATE (icell)
3435 EXIT
3436 END IF
3437 icell(ii, l1, l2, l3) = iat
3438 END DO
3439 IF (ALLOCATED(icell)) EXIT
3440 END DO
3441
3442! duplicate all atoms within boundary layer
3443 laymx = ncx*(2*ll1*ll2 + 2*ll1*ll3 + 2*ll2*ll3 + 4*ll1 + 4*ll2 + 4*ll3 + 8)
3444 nn = nat + laymx
3445 ALLOCATE (rxyz(3, nn), lay(nn))
3446 DO iat = 1, nat
3447 lay(iat) = iat
3448 rxyz(1, iat) = rxyz0(1, iat)
3449 rxyz(2, iat) = rxyz0(2, iat)
3450 rxyz(3, iat) = rxyz0(3, iat)
3451 END DO
3452 il = nat
3453! xy plane
3454 DO l2 = 0, ll2 - 1
3455 DO l1 = 0, ll1 - 1
3456
3457 in = icell(0, l1, l2, 0)
3458 icell(0, l1, l2, ll3) = in
3459 DO ii = 1, in
3460 i = icell(ii, l1, l2, 0)
3461 il = il + 1
3462 IF (il > nn) cpabort("enlarge laymx")
3463 lay(il) = i
3464 icell(ii, l1, l2, ll3) = il
3465 rxyz(1, il) = rxyz(1, i)
3466 rxyz(2, il) = rxyz(2, i)
3467 rxyz(3, il) = rxyz(3, i) + alat(3)
3468 END DO
3469
3470 in = icell(0, l1, l2, ll3 - 1)
3471 icell(0, l1, l2, -1) = in
3472 DO ii = 1, in
3473 i = icell(ii, l1, l2, ll3 - 1)
3474 il = il + 1
3475 IF (il > nn) cpabort("enlarge laymx")
3476 lay(il) = i
3477 icell(ii, l1, l2, -1) = il
3478 rxyz(1, il) = rxyz(1, i)
3479 rxyz(2, il) = rxyz(2, i)
3480 rxyz(3, il) = rxyz(3, i) - alat(3)
3481 END DO
3482
3483 END DO
3484 END DO
3485
3486! yz plane
3487 DO l3 = 0, ll3 - 1
3488 DO l2 = 0, ll2 - 1
3489
3490 in = icell(0, 0, l2, l3)
3491 icell(0, ll1, l2, l3) = in
3492 DO ii = 1, in
3493 i = icell(ii, 0, l2, l3)
3494 il = il + 1
3495 IF (il > nn) cpabort("enlarge laymx")
3496 lay(il) = i
3497 icell(ii, ll1, l2, l3) = il
3498 rxyz(1, il) = rxyz(1, i) + alat(1)
3499 rxyz(2, il) = rxyz(2, i)
3500 rxyz(3, il) = rxyz(3, i)
3501 END DO
3502
3503 in = icell(0, ll1 - 1, l2, l3)
3504 icell(0, -1, l2, l3) = in
3505 DO ii = 1, in
3506 i = icell(ii, ll1 - 1, l2, l3)
3507 il = il + 1
3508 IF (il > nn) cpabort("enlarge laymx")
3509 lay(il) = i
3510 icell(ii, -1, l2, l3) = il
3511 rxyz(1, il) = rxyz(1, i) - alat(1)
3512 rxyz(2, il) = rxyz(2, i)
3513 rxyz(3, il) = rxyz(3, i)
3514 END DO
3515
3516 END DO
3517 END DO
3518
3519! xz plane
3520 DO l3 = 0, ll3 - 1
3521 DO l1 = 0, ll1 - 1
3522
3523 in = icell(0, l1, 0, l3)
3524 icell(0, l1, ll2, l3) = in
3525 DO ii = 1, in
3526 i = icell(ii, l1, 0, l3)
3527 il = il + 1
3528 IF (il > nn) cpabort("enlarge laymx")
3529 lay(il) = i
3530 icell(ii, l1, ll2, l3) = il
3531 rxyz(1, il) = rxyz(1, i)
3532 rxyz(2, il) = rxyz(2, i) + alat(2)
3533 rxyz(3, il) = rxyz(3, i)
3534 END DO
3535
3536 in = icell(0, l1, ll2 - 1, l3)
3537 icell(0, l1, -1, l3) = in
3538 DO ii = 1, in
3539 i = icell(ii, l1, ll2 - 1, l3)
3540 il = il + 1
3541 IF (il > nn) cpabort("enlarge laymx")
3542 lay(il) = i
3543 icell(ii, l1, -1, l3) = il
3544 rxyz(1, il) = rxyz(1, i)
3545 rxyz(2, il) = rxyz(2, i) - alat(2)
3546 rxyz(3, il) = rxyz(3, i)
3547 END DO
3548
3549 END DO
3550 END DO
3551
3552! x axis
3553 DO l1 = 0, ll1 - 1
3554
3555 in = icell(0, l1, 0, 0)
3556 icell(0, l1, ll2, ll3) = in
3557 DO ii = 1, in
3558 i = icell(ii, l1, 0, 0)
3559 il = il + 1
3560 IF (il > nn) cpabort("enlarge laymx")
3561 lay(il) = i
3562 icell(ii, l1, ll2, ll3) = il
3563 rxyz(1, il) = rxyz(1, i)
3564 rxyz(2, il) = rxyz(2, i) + alat(2)
3565 rxyz(3, il) = rxyz(3, i) + alat(3)
3566 END DO
3567
3568 in = icell(0, l1, 0, ll3 - 1)
3569 icell(0, l1, ll2, -1) = in
3570 DO ii = 1, in
3571 i = icell(ii, l1, 0, ll3 - 1)
3572 il = il + 1
3573 IF (il > nn) cpabort("enlarge laymx")
3574 lay(il) = i
3575 icell(ii, l1, ll2, -1) = il
3576 rxyz(1, il) = rxyz(1, i)
3577 rxyz(2, il) = rxyz(2, i) + alat(2)
3578 rxyz(3, il) = rxyz(3, i) - alat(3)
3579 END DO
3580
3581 in = icell(0, l1, ll2 - 1, 0)
3582 icell(0, l1, -1, ll3) = in
3583 DO ii = 1, in
3584 i = icell(ii, l1, ll2 - 1, 0)
3585 il = il + 1
3586 IF (il > nn) cpabort("enlarge laymx")
3587 lay(il) = i
3588 icell(ii, l1, -1, ll3) = il
3589 rxyz(1, il) = rxyz(1, i)
3590 rxyz(2, il) = rxyz(2, i) - alat(2)
3591 rxyz(3, il) = rxyz(3, i) + alat(3)
3592 END DO
3593
3594 in = icell(0, l1, ll2 - 1, ll3 - 1)
3595 icell(0, l1, -1, -1) = in
3596 DO ii = 1, in
3597 i = icell(ii, l1, ll2 - 1, ll3 - 1)
3598 il = il + 1
3599 IF (il > nn) cpabort("enlarge laymx")
3600 lay(il) = i
3601 icell(ii, l1, -1, -1) = il
3602 rxyz(1, il) = rxyz(1, i)
3603 rxyz(2, il) = rxyz(2, i) - alat(2)
3604 rxyz(3, il) = rxyz(3, i) - alat(3)
3605 END DO
3606
3607 END DO
3608
3609! y axis
3610 DO l2 = 0, ll2 - 1
3611
3612 in = icell(0, 0, l2, 0)
3613 icell(0, ll1, l2, ll3) = in
3614 DO ii = 1, in
3615 i = icell(ii, 0, l2, 0)
3616 il = il + 1
3617 IF (il > nn) cpabort("enlarge laymx")
3618 lay(il) = i
3619 icell(ii, ll1, l2, ll3) = il
3620 rxyz(1, il) = rxyz(1, i) + alat(1)
3621 rxyz(2, il) = rxyz(2, i)
3622 rxyz(3, il) = rxyz(3, i) + alat(3)
3623 END DO
3624
3625 in = icell(0, 0, l2, ll3 - 1)
3626 icell(0, ll1, l2, -1) = in
3627 DO ii = 1, in
3628 i = icell(ii, 0, l2, ll3 - 1)
3629 il = il + 1
3630 IF (il > nn) cpabort("enlarge laymx")
3631 lay(il) = i
3632 icell(ii, ll1, l2, -1) = il
3633 rxyz(1, il) = rxyz(1, i) + alat(1)
3634 rxyz(2, il) = rxyz(2, i)
3635 rxyz(3, il) = rxyz(3, i) - alat(3)
3636 END DO
3637
3638 in = icell(0, ll1 - 1, l2, 0)
3639 icell(0, -1, l2, ll3) = in
3640 DO ii = 1, in
3641 i = icell(ii, ll1 - 1, l2, 0)
3642 il = il + 1
3643 IF (il > nn) cpabort("enlarge laymx")
3644 lay(il) = i
3645 icell(ii, -1, l2, ll3) = il
3646 rxyz(1, il) = rxyz(1, i) - alat(1)
3647 rxyz(2, il) = rxyz(2, i)
3648 rxyz(3, il) = rxyz(3, i) + alat(3)
3649 END DO
3650
3651 in = icell(0, ll1 - 1, l2, ll3 - 1)
3652 icell(0, -1, l2, -1) = in
3653 DO ii = 1, in
3654 i = icell(ii, ll1 - 1, l2, ll3 - 1)
3655 il = il + 1
3656 IF (il > nn) cpabort("enlarge laymx")
3657 lay(il) = i
3658 icell(ii, -1, l2, -1) = il
3659 rxyz(1, il) = rxyz(1, i) - alat(1)
3660 rxyz(2, il) = rxyz(2, i)
3661 rxyz(3, il) = rxyz(3, i) - alat(3)
3662 END DO
3663
3664 END DO
3665
3666! z axis
3667 DO l3 = 0, ll3 - 1
3668
3669 in = icell(0, 0, 0, l3)
3670 icell(0, ll1, ll2, l3) = in
3671 DO ii = 1, in
3672 i = icell(ii, 0, 0, l3)
3673 il = il + 1
3674 IF (il > nn) cpabort("enlarge laymx")
3675 lay(il) = i
3676 icell(ii, ll1, ll2, l3) = il
3677 rxyz(1, il) = rxyz(1, i) + alat(1)
3678 rxyz(2, il) = rxyz(2, i) + alat(2)
3679 rxyz(3, il) = rxyz(3, i)
3680 END DO
3681
3682 in = icell(0, ll1 - 1, 0, l3)
3683 icell(0, -1, ll2, l3) = in
3684 DO ii = 1, in
3685 i = icell(ii, ll1 - 1, 0, l3)
3686 il = il + 1
3687 IF (il > nn) cpabort("enlarge laymx")
3688 lay(il) = i
3689 icell(ii, -1, ll2, l3) = il
3690 rxyz(1, il) = rxyz(1, i) - alat(1)
3691 rxyz(2, il) = rxyz(2, i) + alat(2)
3692 rxyz(3, il) = rxyz(3, i)
3693 END DO
3694
3695 in = icell(0, 0, ll2 - 1, l3)
3696 icell(0, ll1, -1, l3) = in
3697 DO ii = 1, in
3698 i = icell(ii, 0, ll2 - 1, l3)
3699 il = il + 1
3700 IF (il > nn) cpabort("enlarge laymx")
3701 lay(il) = i
3702 icell(ii, ll1, -1, l3) = il
3703 rxyz(1, il) = rxyz(1, i) + alat(1)
3704 rxyz(2, il) = rxyz(2, i) - alat(2)
3705 rxyz(3, il) = rxyz(3, i)
3706 END DO
3707
3708 in = icell(0, ll1 - 1, ll2 - 1, l3)
3709 icell(0, -1, -1, l3) = in
3710 DO ii = 1, in
3711 i = icell(ii, ll1 - 1, ll2 - 1, l3)
3712 il = il + 1
3713 IF (il > nn) cpabort("enlarge laymx")
3714 lay(il) = i
3715 icell(ii, -1, -1, l3) = il
3716 rxyz(1, il) = rxyz(1, i) - alat(1)
3717 rxyz(2, il) = rxyz(2, i) - alat(2)
3718 rxyz(3, il) = rxyz(3, i)
3719 END DO
3720
3721 END DO
3722
3723! corners
3724 in = icell(0, 0, 0, 0)
3725 icell(0, ll1, ll2, ll3) = in
3726 DO ii = 1, in
3727 i = icell(ii, 0, 0, 0)
3728 il = il + 1
3729 IF (il > nn) cpabort("enlarge laymx")
3730 lay(il) = i
3731 icell(ii, ll1, ll2, ll3) = il
3732 rxyz(1, il) = rxyz(1, i) + alat(1)
3733 rxyz(2, il) = rxyz(2, i) + alat(2)
3734 rxyz(3, il) = rxyz(3, i) + alat(3)
3735 END DO
3736
3737 in = icell(0, ll1 - 1, 0, 0)
3738 icell(0, -1, ll2, ll3) = in
3739 DO ii = 1, in
3740 i = icell(ii, ll1 - 1, 0, 0)
3741 il = il + 1
3742 IF (il > nn) cpabort("enlarge laymx")
3743 lay(il) = i
3744 icell(ii, -1, ll2, ll3) = il
3745 rxyz(1, il) = rxyz(1, i) - alat(1)
3746 rxyz(2, il) = rxyz(2, i) + alat(2)
3747 rxyz(3, il) = rxyz(3, i) + alat(3)
3748 END DO
3749
3750 in = icell(0, 0, ll2 - 1, 0)
3751 icell(0, ll1, -1, ll3) = in
3752 DO ii = 1, in
3753 i = icell(ii, 0, ll2 - 1, 0)
3754 il = il + 1
3755 IF (il > nn) cpabort("enlarge laymx")
3756 lay(il) = i
3757 icell(ii, ll1, -1, ll3) = il
3758 rxyz(1, il) = rxyz(1, i) + alat(1)
3759 rxyz(2, il) = rxyz(2, i) - alat(2)
3760 rxyz(3, il) = rxyz(3, i) + alat(3)
3761 END DO
3762
3763 in = icell(0, ll1 - 1, ll2 - 1, 0)
3764 icell(0, -1, -1, ll3) = in
3765 DO ii = 1, in
3766 i = icell(ii, ll1 - 1, ll2 - 1, 0)
3767 il = il + 1
3768 IF (il > nn) cpabort("enlarge laymx")
3769 lay(il) = i
3770 icell(ii, -1, -1, ll3) = il
3771 rxyz(1, il) = rxyz(1, i) - alat(1)
3772 rxyz(2, il) = rxyz(2, i) - alat(2)
3773 rxyz(3, il) = rxyz(3, i) + alat(3)
3774 END DO
3775
3776 in = icell(0, 0, 0, ll3 - 1)
3777 icell(0, ll1, ll2, -1) = in
3778 DO ii = 1, in
3779 i = icell(ii, 0, 0, ll3 - 1)
3780 il = il + 1
3781 IF (il > nn) cpabort("enlarge laymx")
3782 lay(il) = i
3783 icell(ii, ll1, ll2, -1) = il
3784 rxyz(1, il) = rxyz(1, i) + alat(1)
3785 rxyz(2, il) = rxyz(2, i) + alat(2)
3786 rxyz(3, il) = rxyz(3, i) - alat(3)
3787 END DO
3788
3789 in = icell(0, ll1 - 1, 0, ll3 - 1)
3790 icell(0, -1, ll2, -1) = in
3791 DO ii = 1, in
3792 i = icell(ii, ll1 - 1, 0, ll3 - 1)
3793 il = il + 1
3794 IF (il > nn) cpabort("enlarge laymx")
3795 lay(il) = i
3796 icell(ii, -1, ll2, -1) = il
3797 rxyz(1, il) = rxyz(1, i) - alat(1)
3798 rxyz(2, il) = rxyz(2, i) + alat(2)
3799 rxyz(3, il) = rxyz(3, i) - alat(3)
3800 END DO
3801
3802 in = icell(0, 0, ll2 - 1, ll3 - 1)
3803 icell(0, ll1, -1, -1) = in
3804 DO ii = 1, in
3805 i = icell(ii, 0, ll2 - 1, ll3 - 1)
3806 il = il + 1
3807 IF (il > nn) cpabort("enlarge laymx")
3808 lay(il) = i
3809 icell(ii, ll1, -1, -1) = il
3810 rxyz(1, il) = rxyz(1, i) + alat(1)
3811 rxyz(2, il) = rxyz(2, i) - alat(2)
3812 rxyz(3, il) = rxyz(3, i) - alat(3)
3813 END DO
3814
3815 in = icell(0, ll1 - 1, ll2 - 1, ll3 - 1)
3816 icell(0, -1, -1, -1) = in
3817 DO ii = 1, in
3818 i = icell(ii, ll1 - 1, ll2 - 1, ll3 - 1)
3819 il = il + 1
3820 IF (il > nn) cpabort("enlarge laymx")
3821 lay(il) = i
3822 icell(ii, -1, -1, -1) = il
3823 rxyz(1, il) = rxyz(1, i) - alat(1)
3824 rxyz(2, il) = rxyz(2, i) - alat(2)
3825 rxyz(3, il) = rxyz(3, i) - alat(3)
3826 END DO
3827
3828 ALLOCATE (lsta(2, nat))
3829 nnbrx = 300
3830 DO
3831 nnbrx = 3*nnbrx/2
3832 ALLOCATE (lstb(nnbrx*nat), rel(5, nnbrx*nat))
3833
3834 indlstx = 0
3835
3836 npr = 1
3837 iam = 0
3838
3839 cut2 = cut**2
3840! assign contiguous portions of the arrays lstb and rel to the threads (this version only contains one thread)
3841 myspace = (nat*nnbrx)/npr
3842 IF (iam == 0) myspaceout = myspace
3843! Verlet list, relative positions
3844 indlst = 0
3845 DO l3 = 0, ll3 - 1
3846 DO l2 = 0, ll2 - 1
3847 DO l1 = 0, ll1 - 1
3848 DO ii = 1, icell(0, l1, l2, l3)
3849 iat = icell(ii, l1, l2, l3)
3850 IF (((iat - 1)*npr)/nat == iam) THEN
3851 lsta(1, iat) = iam*myspace + indlst + 1
3852 CALL sw_sublstiat_l(iat, nn, ncx, ll1, ll2, ll3, l1, l2, l3, myspace, &
3853 rxyz, icell, lstb(iam*myspace + 1), lay, rel(1, iam*myspace + 1), cut2, indlst)
3854 lsta(2, iat) = iam*myspace + indlst
3855 ipb = lsta(1, iat)
3856 ndat = lsta(2, iat) - lsta(1, iat) + 1
3857 END IF
3858 END DO
3859 END DO
3860 END DO
3861 END DO
3862 indlstx = max(indlstx, indlst)
3863
3864 IF (indlstx < myspaceout) EXIT
3865 DEALLOCATE (lstb, rel)
3866 END DO
3867
3868 npr = 1
3869 iam = 0
3870 npjx = 300; npjkx = 6000
3871!end of creating pairlist part------------------------------------------------------------
3872
3873!start energy and force calculation-------------------------------------------------------
3874!set all variables to zero
3875 p = 0.0d0
3876 p3 = 0.0d0
3877 pv3 = 0.0d0
3878 fx(:) = 0.0d0
3879 fy(:) = 0.0d0
3880 fz(:) = 0.0d0
3881 fx3(:) = 0.0d0
3882 fy3(:) = 0.0d0
3883 fz3(:) = 0.0d0
3884!-----------------------------------------------------------------------------------------
3885! triple loop for the 2 and 3-body forces
3886! do 20 i
3887! do 30 j
3888! do 40 k
3889! the pairlists lsta and lstb are used for the perodic boundary conditions
3890!-----------------------------------------------------------------------------------------
3891
3892 DO i = 1, nat
3893 CALL sw_subfeniat_l(i, nat, nnbrx, rel, p, p3, fx, fy, fz, fx3, fy3, fz3, lstb, lsta, isigma, sigma)
3894 END DO
3895!-----------------------------------------------------------------------------------------
3896!*****if necessary,********
3897 DO i = 1, nat
3898 fx(i) = fx(i) + fx3(i)
3899 fy(i) = fy(i) + fy3(i)
3900 fz(i) = fz(i) + fz3(i)
3901 END DO
3902!-----------------------------------------------------------------------------------------
3903
3904 DO i = 1, nat
3905 fxyz(1, i) = fx(i)*esigma
3906 fxyz(2, i) = fy(i)*esigma
3907 fxyz(3, i) = fz(i)*esigma
3908 END DO
3909 etot = (p + p3)*eps
3910 DEALLOCATE (rxyz, icell, lay, lsta, lstb, rel)
3911 END SUBROUTINE eip_stillinger_weber_silicon
3912!End of the force calculation-------------------------------------------------------------
3913
3914! **************************************************************************************************
3915!> \brief ...
3916!> \param c ...
3917!> \return ...
3918! **************************************************************************************************
3919 REAL(8) function f(c)
3920 REAL(8) :: c
3921
3922 REAL(8) :: aa, bb, c4, crainv, ra
3923
3924 parameter(aa=7.049556277d0)
3925 parameter(ra=1.8d0, bb=0.6022245584d0)
3926
3927 IF ((c - ra) < 0.d0) THEN
3928 crainv = 1.0d0/(c - ra)
3929 c4 = c*c*c*c
3930 f = aa*bb*4.0d0/(c4*c)*dexp(crainv) + aa*(bb/(c4) - 1.0d0)*dexp(crainv)*crainv*crainv
3931 ELSE
3932 f = 0.d0
3933 END IF
3934
3935 RETURN
3936 END FUNCTION f
3937
3938! **************************************************************************************************
3939!> \brief ...
3940!> \param d ...
3941!> \return ...
3942! **************************************************************************************************
3943 REAL(8) function pe(d)
3944 REAL(8) :: d
3945
3946 REAL(8), PARAMETER :: aa = 7.049556277d0, bb = 0.6022245584d0, &
3947 ra = 1.8d0
3948
3949 IF ((d - ra) < 0.d0) THEN
3950 pe = aa*(bb/(d*d*d*d) - 1.0d0)*dexp(1.0d0/(d - ra))
3951 ELSE
3952 pe = 0.d0
3953 END IF
3954 RETURN
3955 END FUNCTION pe
3956
3957!------------------------------------------------------------------------------------------
3958! **************************************************************************************************
3959!> \brief ...
3960!> \param iat ...
3961!> \param nn ...
3962!> \param ncx ...
3963!> \param ll1 ...
3964!> \param ll2 ...
3965!> \param ll3 ...
3966!> \param l1 ...
3967!> \param l2 ...
3968!> \param l3 ...
3969!> \param myspace ...
3970!> \param rxyz ...
3971!> \param icell ...
3972!> \param lstb ...
3973!> \param lay ...
3974!> \param rel ...
3975!> \param cut2 ...
3976!> \param indlst ...
3977! **************************************************************************************************
3978 SUBROUTINE sw_sublstiat_l(iat, nn, ncx, ll1, ll2, ll3, l1, l2, l3, myspace, &
3979 rxyz, icell, lstb, lay, rel, cut2, indlst)
3980! finds the neighbours of atom iat (specified by lsta and lstb) and and
3981! the relative position rel of iat with respect to these neighbours
3982 INTEGER :: iat, nn, ncx, ll1, ll2, ll3, l1, l2, l3, &
3983 myspace
3984 REAL(8) :: rxyz(3, nn)
3985 INTEGER :: icell(0:ncx, -1:ll1, -1:ll2, -1:ll3), lstb(0:myspace - 1), lay(nn)
3986 REAL(8) :: rel(5, 0:myspace - 1), cut2
3987 INTEGER :: indlst
3988
3989 INTEGER :: jat, jj, k1, k2, k3
3990 REAL(8) :: rr2, tt, tti, xrel, yrel, zrel
3991
3992 DO k3 = l3 - 1, l3 + 1
3993 DO k2 = l2 - 1, l2 + 1
3994 DO k1 = l1 - 1, l1 + 1
3995 DO jj = 1, icell(0, k1, k2, k3)
3996 jat = icell(jj, k1, k2, k3)
3997 IF (jat == iat) cycle
3998 xrel = rxyz(1, iat) - rxyz(1, jat)
3999 yrel = rxyz(2, iat) - rxyz(2, jat)
4000 zrel = rxyz(3, iat) - rxyz(3, jat)
4001 rr2 = xrel**2 + yrel**2 + zrel**2
4002 IF (rr2 <= cut2) THEN
4003 indlst = min(indlst, myspace - 1)
4004 lstb(indlst) = lay(jat)
4005! write(6,*) 'iat,indlst,lay(jat)',iat,indlst,lay(jat)
4006 tt = sqrt(rr2)
4007 tti = 1.d0/tt
4008 rel(1, indlst) = xrel*tti
4009 rel(2, indlst) = yrel*tti
4010 rel(3, indlst) = zrel*tti
4011 rel(4, indlst) = tt
4012 rel(5, indlst) = tti
4013 indlst = indlst + 1
4014 END IF
4015 END DO
4016 END DO
4017 END DO
4018 END DO
4019
4020 RETURN
4021 END SUBROUTINE sw_sublstiat_l
4022
4023! **************************************************************************************************
4024!> \brief ...
4025!> \param i ...
4026!> \param nat ...
4027!> \param nnbrx ...
4028!> \param rel ...
4029!> \param p ...
4030!> \param p3 ...
4031!> \param fx ...
4032!> \param fy ...
4033!> \param fz ...
4034!> \param fx3 ...
4035!> \param fy3 ...
4036!> \param fz3 ...
4037!> \param lstb ...
4038!> \param lsta ...
4039!> \param isigma ...
4040!> \param sigma ...
4041! **************************************************************************************************
4042 SUBROUTINE sw_subfeniat_l(i, nat, nnbrx, rel, p, p3, fx, fy, fz, fx3, fy3, fz3, lstb, lsta, isigma, sigma)
4043 INTEGER, INTENT(IN) :: i, nat, nnbrx
4044 REAL(8), INTENT(IN) :: rel(5, nnbrx*nat)
4045 REAL(8), INTENT(INOUT) :: p, p3, fx(nat), fy(nat), fz(nat), &
4046 fx3(nat), fy3(nat), fz3(nat)
4047 INTEGER, INTENT(IN) :: lstb(nnbrx*nat), lsta(2, nat)
4048 REAL(8), INTENT(IN) :: isigma, sigma
4049
4050 INTEGER :: ipb, ipe, j, k, l, m, nij
4051 REAL(8) :: aa, bb, c4, cosijk, cosijk3, cosikj, cosikj3, cosjik, cosjik3, crainv, force, &
4052 gam, hi, hixij, hixij0, hixij1, hixik, hixik0, hixik1, hiyij, hiyij0, hiyij1, hiyik, &
4053 hiyik0, hiyik1, hizij, hizij0, hizij1, hizik, hizik0, hizik1, hj, hjxij, hjxij0, hjxij1, &
4054 hjxjk, hjxjk0, hjxjk1, hjyij, hjyij0, hjyij1, hjyjk, hjyjk0, hjyjk1, hjzij, hjzij0, &
4055 hjzij1, hjzjk, hjzjk0, hjzjk1, hk, hkxik, hkxik0, hkxik1, hkxkj, hkxkj0, hkxkj1, hkyik, &
4056 hkyik0, hkyik1, hkykj, hkykj0, hkykj1, hkzik, hkzik0, hkzik1, hkzkj, hkzkj0, hkzkj1, &
4057 invrij, invrija, invrik, invrika, invrjk, invrjka, ra, ramda, refi, refj
4058 REAL(8) :: refk, rij, rija, rik, rika, rjk, rjka, xij, xik, xjk, yij, yik, yjk, zij, zik, zjk
4059
4060 parameter(gam=1.2d0, ramda=21.0d0, aa=7.049556277d0)
4061 parameter(ra=1.8d0, bb=0.6022245584d0)
4062
4063 ipb = lsta(1, i)
4064 ipe = lsta(2, i)
4065
4066 DO l = ipb, ipe
4067 j = lstb(l)
4068 IF (j <= i) cycle
4069 nij = 0
4070 rij = rel(4, l)*isigma
4071 invrij = rel(5, l)*sigma
4072 xij = rel(1, l)*rij
4073 yij = rel(2, l)*rij
4074 zij = rel(3, l)*rij
4075
4076 IF (rij >= 2.d0*ra) cycle
4077 IF (rij < ra) THEN
4078 crainv = 1.0d0/(rij - ra)
4079 c4 = rij*rij*rij*rij
4080 force = aa*bb*4.0d0/(c4*rij)*dexp(crainv) + aa*(bb/(c4) - 1.0d0)*dexp(crainv)*crainv*crainv
4081
4082 fx(i) = force*xij*invrij + fx(i)
4083 fy(i) = force*yij*invrij + fy(i)
4084 fz(i) = force*zij*invrij + fz(i)
4085
4086 fx(j) = -force*xij*invrij + fx(j)
4087 fy(j) = -force*yij*invrij + fy(j)
4088 fz(j) = -force*zij*invrij + fz(j)
4089
4090 p = p + aa*(bb/(rij*rij*rij*rij) - 1.0d0)*dexp(1.0d0/(rij - ra))
4091
4092 nij = 1
4093 END IF
4094
4095 DO m = ipb, ipe
4096 k = lstb(m)
4097 IF (k <= j) cycle
4098 invrik = rel(5, m)*sigma
4099 rik = rel(4, m)*isigma
4100 xik = rel(1, m)*rik
4101 yik = rel(2, m)*rik
4102 zik = rel(3, m)*rik
4103
4104 IF ((rik >= ra) .AND. (nij == 0)) cycle
4105
4106 xjk = xik - xij
4107 yjk = yik - yij
4108 zjk = zik - zij
4109
4110 rjk = sqrt(xjk*xjk + yjk*yjk + zjk*zjk)
4111 invrjk = 1.d0/rjk
4112
4113 IF ((rjk >= ra) .AND. (nij == 0)) cycle
4114 cosjik = (xij*xik + yij*yik + zij*zik)*(invrij*invrik)
4115 cosijk = (-xij*xjk - yij*yjk - zij*zjk)*(invrij*invrjk)
4116 cosikj = (xik*xjk + yik*yjk + zik*zjk)*(invrik*invrjk)
4117 cosjik3 = cosjik + 1.0d0/3.0d0
4118 cosijk3 = cosijk + 1.0d0/3.0d0
4119 cosikj3 = cosikj + 1.0d0/3.0d0
4120
4121 rija = rij - ra
4122 rika = rik - ra
4123 rjka = rjk - ra
4124
4125 invrija = 1.d0/rija
4126 invrika = 1.d0/rika
4127 invrjka = 1.d0/rjka
4128
4129 IF (rija >= 0.0d0) THEN
4130 refi = 0.0d0
4131 refj = 0.0d0
4132 refk = ramda*exp(gam*invrika + gam*invrjka)
4133 ELSE IF ((rija < 0.0d0) .AND. (rika < 0.0d0)) THEN
4134 IF (rjka < 0.0d0) THEN
4135 refi = ramda*exp(gam*invrija + gam*invrika)
4136 refj = ramda*exp(gam*invrija + gam*invrjka)
4137 refk = ramda*exp(gam*invrika + gam*invrjka)
4138 ELSE
4139 refi = ramda*exp(gam*invrija + gam*invrika)
4140 refj = 0.0d0
4141 refk = 0.0d0
4142 END IF
4143 ELSE IF ((rija < 0.0d0) .AND. (rjka < 0.0d0)) THEN
4144 refi = 0.0d0
4145 refj = ramda*exp(gam*invrija + gam*invrjka)
4146 refk = 0.0d0
4147 ELSE
4148 cycle
4149 END IF
4150
4151 hi = refi*cosjik3*cosjik3
4152 hj = refj*cosijk3*cosijk3
4153 hk = refk*cosikj3*cosikj3
4154 p3 = p3 + hi + hj + hk
4155
4156 hixij0 = 2.0d0*(xik*invrik - xij*cosjik*invrij)
4157 hixij1 = gam*xij*cosjik3*(invrija*invrija)
4158 hixij = refi*cosjik3*(hixij0 - hixij1)*invrij
4159 hixik0 = 2.0d0*(xij*invrij - xik*cosjik*invrik)
4160 hixik1 = gam*xik*cosjik3*(invrika*invrika)
4161 hixik = refi*cosjik3*(hixik0 - hixik1)*invrik
4162 hjxij0 = 2.0d0*(-xjk*invrjk - xij*cosijk*invrij)
4163 hjxij1 = gam*xij*cosijk3*(invrija*invrija)
4164 hjxij = refj*cosijk3*(hjxij0 - hjxij1)*invrij
4165 hkxik0 = 2.0d0*(xjk*invrjk - xik*cosikj*invrik)
4166 hkxik1 = gam*xik*cosikj3*(invrika*invrika)
4167 hkxik = refk*cosikj3*(hkxik0 - hkxik1)*invrik
4168 hjxjk0 = 2.0d0*(-xij*invrij - xjk*cosijk*invrjk)
4169 hjxjk1 = gam*xjk*cosijk3*(invrjka*invrjka)
4170 hjxjk = refj*cosijk3*(hjxjk0 - hjxjk1)*invrjk
4171 hkxkj0 = 2.0d0*(-xik*invrik + xjk*cosikj*invrjk)
4172 hkxkj1 = gam*xjk*cosikj3*(invrjka*invrjka)
4173 hkxkj = refk*cosikj3*(hkxkj0 + hkxkj1)*invrjk
4174
4175 hiyij0 = 2.0d0*(yik*invrik - yij*cosjik*invrij)
4176 hiyij1 = gam*yij*cosjik3*(invrija*invrija)
4177 hiyij = refi*cosjik3*(hiyij0 - hiyij1)*invrij
4178 hiyik0 = 2.0d0*(yij*invrij - yik*cosjik*invrik)
4179 hiyik1 = gam*yik*cosjik3*(invrika*invrika)
4180 hiyik = refi*cosjik3*(hiyik0 - hiyik1)*invrik
4181 hjyij0 = 2.0d0*(-yjk*invrjk - yij*cosijk*invrij)
4182 hjyij1 = gam*yij*cosijk3*(invrija*invrija)
4183 hjyij = refj*cosijk3*(hjyij0 - hjyij1)*invrij
4184 hkyik0 = 2.0d0*(yjk*invrjk - yik*cosikj*invrik)
4185 hkyik1 = gam*yik*cosikj3*(invrika*invrika)
4186 hkyik = refk*cosikj3*(hkyik0 - hkyik1)*invrik
4187 hjyjk0 = 2.0d0*(-yij*invrij - yjk*cosijk*invrjk)
4188 hjyjk1 = gam*yjk*cosijk3*(invrjka*invrjka)
4189 hjyjk = refj*cosijk3*(hjyjk0 - hjyjk1)*invrjk
4190 hkykj0 = 2.0d0*(-yik*invrik + yjk*cosikj*invrjk)
4191 hkykj1 = gam*yjk*cosikj3*(invrjka*invrjka)
4192 hkykj = refk*cosikj3*(hkykj0 + hkykj1)*invrjk
4193
4194 hizij0 = 2.0d0*(zik*invrik - zij*cosjik*invrij)
4195 hizij1 = gam*zij*cosjik3*(invrija*invrija)
4196 hizij = refi*cosjik3*(hizij0 - hizij1)*invrij
4197 hizik0 = 2.0d0*(zij*invrij - zik*cosjik*invrik)
4198 hizik1 = gam*zik*cosjik3*(invrika*invrika)
4199 hizik = refi*cosjik3*(hizik0 - hizik1)*invrik
4200 hjzij0 = 2.0d0*(-zjk*invrjk - zij*cosijk*invrij)
4201 hjzij1 = gam*zij*cosijk3*(invrija*invrija)
4202 hjzij = refj*cosijk3*(hjzij0 - hjzij1)*invrij
4203 hkzik0 = 2.0d0*(zjk*invrjk - zik*cosikj*invrik)
4204 hkzik1 = gam*zik*cosikj3*(invrika*invrika)
4205 hkzik = refk*cosikj3*(hkzik0 - hkzik1)*invrik
4206 hjzjk0 = 2.0d0*(-zij*invrij - zjk*cosijk*invrjk)
4207 hjzjk1 = gam*zjk*cosijk3*(invrjka*invrjka)
4208 hjzjk = refj*cosijk3*(hjzjk0 - hjzjk1)*invrjk
4209 hkzkj0 = 2.0d0*(-zik*invrik + zjk*cosikj*invrjk)
4210 hkzkj1 = gam*zjk*cosikj3*(invrjka*invrjka)
4211 hkzkj = refk*cosikj3*(hkzkj0 + hkzkj1)*invrjk
4212
4213 fx3(i) = fx3(i) - hixij - hixik - hjxij - hkxik
4214 fy3(i) = fy3(i) - hiyij - hiyik - hjyij - hkyik
4215 fz3(i) = fz3(i) - hizij - hizik - hjzij - hkzik
4216
4217 fx3(j) = fx3(j) + hixij + hjxij - hjxjk + hkxkj
4218 fy3(j) = fy3(j) + hiyij + hjyij - hjyjk + hkykj
4219 fz3(j) = fz3(j) + hizij + hjzij - hjzjk + hkzkj
4220
4221 fx3(k) = fx3(k) + hixik + hkxik - hkxkj + hjxjk
4222 fy3(k) = fy3(k) + hiyik + hkyik - hkykj + hjyjk
4223 fz3(k) = fz3(k) + hizik + hkzik - hkzkj + hjzjk
4224 END DO
4225 END DO
4226 END SUBROUTINE sw_subfeniat_l
4227
4228! **************************************************************************************************
4229!> \brief ...
4230!> \param nat ...
4231!> \param alat ...
4232!> \param rxyz ...
4233!> \param fxyz ...
4234!> \param etot ...
4235!> \param count ...
4236! **************************************************************************************************
4237 SUBROUTINE eip_tersoff_silicon(nat, alat, rxyz, fxyz, etot, count)
4238!*****************************************************************************************
4239! This subroutine evaluates the Tersoff Silicon potential with linear scaling
4240! COPYRIGHT
4241! Copyright (C) 2009 AIST, UNIBAS
4242! This file is distributed under the terms of the
4243! GNU General Public License, see
4244! http://www.gnu.org/copyleft/gpl.txt .
4245!
4246! Implementation: Original version was written by Kengo Nishio, AIST Tsukuba (JP)
4247! Improved by M. Amsler, S. Goedecker, Basel University (CH), 2009
4248!
4249! Note:
4250! Parameters and functional form from PRL 61, 2879 (1988) and PRB 39, 5566 (1989)
4251!
4252! Input:
4253! nat, integer: the number of atoms
4254! alat, real(8), dim(3) : the three edges of the orthoromic simulation cell, periodic boundaries are applied
4255! and atoms outside the cell will be brought back into the box
4256! rxyz, real(8), dim(3,nat) : the xyz cartesian components of the atomic positions in Angstroem
4257!
4258! Output:
4259! fxyz, real(8), dim(3,nat): the xyz cartesian forces in on the corresponding atomic components n eV/A
4260! etot, real(8) : total potential energy, 2-body and 3-body, in eV
4261! count, real(8) : increased by 1.d0 at each call of this subroutine,
4262! needs to be initialized to 0.d0 before calling this routine for the first time
4263!*****************************************************************************************
4264 INTEGER :: nat
4265 REAL(8) :: alat(3), rxyz(3, nat), fxyz(3, nat), &
4266 etot, count
4267
4268 INTEGER :: iat, nnmax, npmax
4269 INTEGER, ALLOCATABLE, DIMENSION(:) :: kinds, lstb
4270 INTEGER, ALLOCATABLE, DIMENSION(:, :) :: lsta
4271 REAL(8) :: uatot, urtot, xbox, ybox, zbox
4272 REAL(8), ALLOCATABLE, DIMENSION(:) :: dkeij, uadurdf, xyzrrefdf
4273 REAL(8), DIMENSION(1:2) :: bcsq, co_bcd, dsq, h, pmass, pn
4274 REAL(8), DIMENSION(1:2, 1:2) :: ala, alr, ca, cr, r1, r2, x
4275
4276 INTEGER:: nnbrx, nnbrxt
4277 INTEGER :: i
4278
4279 count = count + 1.d0
4280
4281 DO iat = 1, nat
4282 rxyz(1, iat) = modulo(modulo(rxyz(1, iat), alat(1)), alat(1))
4283 rxyz(2, iat) = modulo(modulo(rxyz(2, iat), alat(2)), alat(2))
4284 rxyz(3, iat) = modulo(modulo(rxyz(3, iat), alat(3)), alat(3))
4285 END DO
4286
4287 ALLOCATE (kinds(1:nat))
4288
4289 nnbrx = 24
4290 nnbrxt = 3*nnbrx/2
4291 nnmax = nnbrxt*nat
4292 npmax = nnbrxt*nat
4293 ALLOCATE (lsta(2, nat), lstb(nnbrxt*nat))
4294 ALLOCATE (xyzrrefdf(1:6*npmax), uadurdf(1:3*npmax), dkeij(1:3*nnmax))
4295
4296 DO i = 1, nat
4297 kinds(i) = 2 !Since all atoms are Si, all of kind 2
4298 END DO
4299 fxyz = 0.0d0
4300 xbox = alat(1); ybox = alat(2); zbox = alat(3)
4301 CALL tersoff_parameters(r1, r2, cr, ca, alr, ala, x, pn, co_bcd, bcsq, dsq, h, pmass)
4302 CALL tersoff_pairlist_energy_forces(nat, npmax, nnmax, xbox, ybox, zbox, kinds, rxyz, r1, r2, cr, &
4303 ca, alr, ala, x, xyzrrefdf, uadurdf, urtot, lsta, lstb, nnbrx, &
4304 pn, co_bcd, bcsq, dsq, h, fxyz, uatot, dkeij)
4305 etot = urtot + uatot
4306 DEALLOCATE (kinds, xyzrrefdf, uadurdf, dkeij, lsta, lstb)
4307 END SUBROUTINE eip_tersoff_silicon
4308!-----------------------------------------------------------------------------------------
4309! **************************************************************************************************
4310!> \brief ...
4311!> \param R1 ...
4312!> \param R2 ...
4313!> \param Cr ...
4314!> \param Ca ...
4315!> \param alr ...
4316!> \param ala ...
4317!> \param X ...
4318!> \param Pn ...
4319!> \param Co_bcd ...
4320!> \param bcsq ...
4321!> \param dsq ...
4322!> \param h ...
4323!> \param Pmass ...
4324! **************************************************************************************************
4325 SUBROUTINE tersoff_parameters(R1, R2, Cr, Ca, alr, ala, X, Pn, Co_bcd, bcsq, dsq, h, Pmass)
4326
4327 REAL(8), DIMENSION(1:2, 1:2), INTENT(out) :: r1, r2, cr, ca, alr, ala, x
4328 REAL(8), DIMENSION(1:2), INTENT(out) :: pn, co_bcd, bcsq, dsq, h, pmass
4329
4330 REAL(8), PARAMETER :: c_ala = 2.2119d0, c_alr = 3.4879d0, c_b = 1.5724d-7, c_c = 3.8049d4, &
4331 c_ca = 3.4674d2, c_cr = 1.3936d3, c_d = 4.3484d0, c_h = -5.7058d-1, c_mass = 12.0d0, &
4332 c_n = 7.2751d-1, c_r1 = 1.8d0, c_r2 = 2.1d0, si_ala = 1.7322d0, si_alr = 2.4799d0, &
4333 si_b = 1.1000d-6, si_c = 1.0039d5, si_ca = 4.7118d2, si_cr = 1.8308d3, si_d = 1.6217d1, &
4334 si_h = -5.9825d-1, si_mass = 28.0855d0, si_n = 7.8734d-1, si_r1 = 2.7d0, si_r2 = 3.3d0
4335
4336!Parameter for carbon, not used in this version
4337!Parameter for carbon, not used in this version
4338!Parameter for carbon, not used in this version
4339!Parameter for carbon, not used in this version
4340!Parameter for carbon, not used in this version
4341!Parameter for carbon, not used in this version
4342!Parameter for carbon, not used in this version
4343!Parameter for carbon, not used in this version
4344!Parameter for carbon, not used in this version
4345!Parameter for carbon, not used in this version
4346!Parameter for carbon, not used in this version
4347!Parameter for carbon, not used in this version
4348!Increased Cutoff, originally 3.0d0
4349
4350 cr(1, 1) = c_cr
4351 cr(2, 2) = si_cr
4352 cr(1, 2) = dsqrt(cr(1, 1)*cr(2, 2))
4353 cr(2, 1) = cr(1, 2)
4354
4355 ca(1, 1) = c_ca
4356 ca(2, 2) = si_ca
4357 ca(1, 2) = dsqrt(ca(1, 1)*ca(2, 2))
4358 ca(2, 1) = ca(1, 2)
4359
4360 r1(1, 1) = c_r1
4361 r1(2, 2) = si_r1
4362 r1(1, 2) = dsqrt(r1(1, 1)*r1(2, 2))
4363 r1(2, 1) = r1(1, 2)
4364
4365 r2(1, 1) = c_r2
4366 r2(2, 2) = si_r2
4367 r2(1, 2) = dsqrt(r2(1, 1)*r2(2, 2))
4368 r2(2, 1) = r2(1, 2)
4369
4370 x(1, 1) = 1.0d0
4371 x(2, 2) = 1.0d0
4372 x(1, 2) = 0.9776d0
4373 x(2, 1) = 0.9776d0
4374
4375 alr(1, 1) = c_alr
4376 alr(2, 2) = si_alr
4377 alr(1, 2) = 0.5d0*(alr(1, 1) + alr(2, 2))
4378 alr(2, 1) = alr(1, 2)
4379
4380 ala(1, 1) = c_ala
4381 ala(2, 2) = si_ala
4382 ala(1, 2) = 0.5d0*(ala(1, 1) + ala(2, 2))
4383 ala(2, 1) = ala(1, 2)
4384
4385 pn(1) = c_n
4386 pn(2) = si_n
4387
4388 co_bcd(1) = c_b*(1.0d0 + c_c*c_c/(c_d*c_d))
4389 co_bcd(2) = si_b*(1.0d0 + si_c*si_c/(si_d*si_d))
4390
4391 bcsq(1) = c_b*c_c*c_c
4392 bcsq(2) = si_b*si_c*si_c
4393
4394 dsq(1) = c_d*c_d
4395 dsq(2) = si_d*si_d
4396
4397 h(1) = c_h
4398 h(2) = si_h
4399
4400 pmass(1) = c_mass
4401 pmass(2) = si_mass
4402
4403 RETURN
4404 END SUBROUTINE tersoff_parameters
4405!-----------------------------------------------------------------------------------------
4406! **************************************************************************************************
4407!> \brief ...
4408!> \param Nmol ...
4409!> \param Npmax ...
4410!> \param NNmax ...
4411!> \param xbox ...
4412!> \param ybox ...
4413!> \param zbox ...
4414!> \param Kinds ...
4415!> \param R ...
4416!> \param R1 ...
4417!> \param R2 ...
4418!> \param Cr ...
4419!> \param Ca ...
4420!> \param alr ...
4421!> \param ala ...
4422!> \param X ...
4423!> \param XYZRrefdf ...
4424!> \param UadUrdf ...
4425!> \param Urtot ...
4426!> \param lsta ...
4427!> \param lstb ...
4428!> \param nnbrx ...
4429!> \param Pn ...
4430!> \param Co_bcd ...
4431!> \param bcsq ...
4432!> \param dsq ...
4433!> \param h ...
4434!> \param F ...
4435!> \param Uatot ...
4436!> \param dkEij ...
4437! **************************************************************************************************
4438 SUBROUTINE tersoff_pairlist_energy_forces(Nmol, Npmax, NNmax, xbox, ybox, zbox, Kinds, R, R1, R2, Cr, Ca, alr, ala, X, &
4439 XYZRrefdf, UadUrdf, Urtot, lsta, lstb, nnbrx, &
4440 Pn, Co_bcd, bcsq, dsq, h, F, Uatot, dkEij)
4441 INTEGER, INTENT(in) :: nmol, npmax, nnmax
4442 REAL(8), INTENT(in) :: xbox, ybox, zbox
4443 INTEGER, DIMENSION(1:Nmol), INTENT(in) :: kinds
4444 REAL(8), DIMENSION(1:3*Nmol), INTENT(in) :: r
4445 REAL(8), DIMENSION(1:2, 1:2), INTENT(in) :: r1, r2, cr, ca, alr, ala, x
4446 REAL(8), DIMENSION(1:6*Npmax), INTENT(out) :: xyzrrefdf
4447 REAL(8), DIMENSION(1:3*Npmax), INTENT(out) :: uadurdf
4448 REAL(8), INTENT(out) :: urtot
4449 INTEGER :: lsta(2, nmol)
4450 INTEGER, INTENT(inout) :: nnbrx
4451 INTEGER :: lstb(nnbrx*nmol)
4452 REAL(8), DIMENSION(1:2), INTENT(in) :: pn, co_bcd, bcsq, dsq, h
4453 REAL(8), DIMENSION(1:3*Nmol), INTENT(out) :: f
4454 REAL(8), INTENT(out) :: uatot
4455 REAL(8), DIMENSION(1:3*NNmax) :: dkeij
4456
4457 INTEGER :: i, iam, iat, ii, il, in, indlst, indlstx, ipb, istopg, jat, l1, l2, l3, laymx, &
4458 ll1, ll2, ll3, myspace, myspaceout, nat, ncx, ndat, nn, npjkx, npjx, npr, nptot
4459 INTEGER, ALLOCATABLE, DIMENSION(:) :: lay
4460 INTEGER, ALLOCATABLE, DIMENSION(:, :, :, :) :: icell
4461 REAL(8) :: alat(3), cut, cut2, pi, rlc1i, rlc2i, &
4462 rlc3i, rxyz0(3, nmol), xhalf, yhalf, &
4463 zhalf
4464 REAL(8), ALLOCATABLE, DIMENSION(:, :) :: rel, rxyz
4465
4466 pi = dacos(-1.0d0)
4467
4468 xhalf = 0.5d0*xbox
4469 yhalf = 0.5d0*ybox
4470 zhalf = 0.5d0*zbox
4471
4472 nat = nmol
4473
4474 alat(1) = xbox
4475 alat(2) = ybox
4476 alat(3) = zbox
4477
4478 DO iat = 1, nat
4479 jat = 3*(iat - 1)
4480 rxyz0(1, iat) = r(jat + 1)
4481 rxyz0(2, iat) = r(jat + 2)
4482 rxyz0(3, iat) = r(jat + 3)
4483 END DO
4484
4485 cut = r2(2, 2) - 1.d-9
4486
4487! linear scaling calculation of verlet list
4488 ll1 = int(alat(1)/cut)
4489 IF (ll1 < 1) cpabort("alat(1) too small")
4490 ll2 = int(alat(2)/cut)
4491 IF (ll2 < 1) cpabort("alat(2) too small")
4492 ll3 = int(alat(3)/cut)
4493 IF (ll3 < 1) cpabort("alat(3) too small")
4494
4495! determine number of threadsi (this version is only singlethreaded)
4496 npr = 1
4497! linear scaling calculation of verlet list
4498
4499 ncx = 8
4500 DO
4501 ncx = ncx*2
4502 ALLOCATE (icell(0:ncx, -1:ll1, -1:ll2, -1:ll3))
4503 icell(0, :, :, :) = 0
4504 rlc1i = ll1/alat(1)
4505 rlc2i = ll2/alat(2)
4506 rlc3i = ll3/alat(3)
4507
4508 DO iat = 1, nat
4509 l1 = int(rxyz0(1, iat)*rlc1i)
4510 l2 = int(rxyz0(2, iat)*rlc2i)
4511 l3 = int(rxyz0(3, iat)*rlc3i)
4512
4513 ii = icell(0, l1, l2, l3)
4514 ii = ii + 1
4515 icell(0, l1, l2, l3) = ii
4516 IF (ii > ncx) THEN
4517 DEALLOCATE (icell)
4518 EXIT
4519 END IF
4520 icell(ii, l1, l2, l3) = iat
4521 END DO
4522 IF (ALLOCATED(icell)) EXIT
4523 END DO
4524
4525! duplicate all atoms within boundary layer
4526 laymx = ncx*(2*ll1*ll2 + 2*ll1*ll3 + 2*ll2*ll3 + 4*ll1 + 4*ll2 + 4*ll3 + 8)
4527 nn = nat + laymx
4528 ALLOCATE (rxyz(3, nn), lay(nn))
4529 DO iat = 1, nat
4530 lay(iat) = iat
4531 rxyz(1, iat) = rxyz0(1, iat)
4532 rxyz(2, iat) = rxyz0(2, iat)
4533 rxyz(3, iat) = rxyz0(3, iat)
4534 END DO
4535 il = nat
4536! xy plane
4537 DO l2 = 0, ll2 - 1
4538 DO l1 = 0, ll1 - 1
4539
4540 in = icell(0, l1, l2, 0)
4541 icell(0, l1, l2, ll3) = in
4542 DO ii = 1, in
4543 i = icell(ii, l1, l2, 0)
4544 il = il + 1
4545 IF (il > nn) cpabort("enlarge laymx")
4546 lay(il) = i
4547 icell(ii, l1, l2, ll3) = il
4548 rxyz(1, il) = rxyz(1, i)
4549 rxyz(2, il) = rxyz(2, i)
4550 rxyz(3, il) = rxyz(3, i) + alat(3)
4551 END DO
4552
4553 in = icell(0, l1, l2, ll3 - 1)
4554 icell(0, l1, l2, -1) = in
4555 DO ii = 1, in
4556 i = icell(ii, l1, l2, ll3 - 1)
4557 il = il + 1
4558 IF (il > nn) cpabort("enlarge laymx")
4559 lay(il) = i
4560 icell(ii, l1, l2, -1) = il
4561 rxyz(1, il) = rxyz(1, i)
4562 rxyz(2, il) = rxyz(2, i)
4563 rxyz(3, il) = rxyz(3, i) - alat(3)
4564 END DO
4565
4566 END DO
4567 END DO
4568
4569! yz plane
4570 DO l3 = 0, ll3 - 1
4571 DO l2 = 0, ll2 - 1
4572
4573 in = icell(0, 0, l2, l3)
4574 icell(0, ll1, l2, l3) = in
4575 DO ii = 1, in
4576 i = icell(ii, 0, l2, l3)
4577 il = il + 1
4578 IF (il > nn) cpabort("enlarge laymx")
4579 lay(il) = i
4580 icell(ii, ll1, l2, l3) = il
4581 rxyz(1, il) = rxyz(1, i) + alat(1)
4582 rxyz(2, il) = rxyz(2, i)
4583 rxyz(3, il) = rxyz(3, i)
4584 END DO
4585
4586 in = icell(0, ll1 - 1, l2, l3)
4587 icell(0, -1, l2, l3) = in
4588 DO ii = 1, in
4589 i = icell(ii, ll1 - 1, l2, l3)
4590 il = il + 1
4591 IF (il > nn) cpabort("enlarge laymx")
4592 lay(il) = i
4593 icell(ii, -1, l2, l3) = il
4594 rxyz(1, il) = rxyz(1, i) - alat(1)
4595 rxyz(2, il) = rxyz(2, i)
4596 rxyz(3, il) = rxyz(3, i)
4597 END DO
4598
4599 END DO
4600 END DO
4601
4602! xz plane
4603 DO l3 = 0, ll3 - 1
4604 DO l1 = 0, ll1 - 1
4605
4606 in = icell(0, l1, 0, l3)
4607 icell(0, l1, ll2, l3) = in
4608 DO ii = 1, in
4609 i = icell(ii, l1, 0, l3)
4610 il = il + 1
4611 IF (il > nn) cpabort("enlarge laymx")
4612 lay(il) = i
4613 icell(ii, l1, ll2, l3) = il
4614 rxyz(1, il) = rxyz(1, i)
4615 rxyz(2, il) = rxyz(2, i) + alat(2)
4616 rxyz(3, il) = rxyz(3, i)
4617 END DO
4618
4619 in = icell(0, l1, ll2 - 1, l3)
4620 icell(0, l1, -1, l3) = in
4621 DO ii = 1, in
4622 i = icell(ii, l1, ll2 - 1, l3)
4623 il = il + 1
4624 IF (il > nn) cpabort("enlarge laymx")
4625 lay(il) = i
4626 icell(ii, l1, -1, l3) = il
4627 rxyz(1, il) = rxyz(1, i)
4628 rxyz(2, il) = rxyz(2, i) - alat(2)
4629 rxyz(3, il) = rxyz(3, i)
4630 END DO
4631
4632 END DO
4633 END DO
4634
4635! x axis
4636 DO l1 = 0, ll1 - 1
4637
4638 in = icell(0, l1, 0, 0)
4639 icell(0, l1, ll2, ll3) = in
4640 DO ii = 1, in
4641 i = icell(ii, l1, 0, 0)
4642 il = il + 1
4643 IF (il > nn) cpabort("enlarge laymx")
4644 lay(il) = i
4645 icell(ii, l1, ll2, ll3) = il
4646 rxyz(1, il) = rxyz(1, i)
4647 rxyz(2, il) = rxyz(2, i) + alat(2)
4648 rxyz(3, il) = rxyz(3, i) + alat(3)
4649 END DO
4650
4651 in = icell(0, l1, 0, ll3 - 1)
4652 icell(0, l1, ll2, -1) = in
4653 DO ii = 1, in
4654 i = icell(ii, l1, 0, ll3 - 1)
4655 il = il + 1
4656 IF (il > nn) cpabort("enlarge laymx")
4657 lay(il) = i
4658 icell(ii, l1, ll2, -1) = il
4659 rxyz(1, il) = rxyz(1, i)
4660 rxyz(2, il) = rxyz(2, i) + alat(2)
4661 rxyz(3, il) = rxyz(3, i) - alat(3)
4662 END DO
4663
4664 in = icell(0, l1, ll2 - 1, 0)
4665 icell(0, l1, -1, ll3) = in
4666 DO ii = 1, in
4667 i = icell(ii, l1, ll2 - 1, 0)
4668 il = il + 1
4669 IF (il > nn) cpabort("enlarge laymx")
4670 lay(il) = i
4671 icell(ii, l1, -1, ll3) = il
4672 rxyz(1, il) = rxyz(1, i)
4673 rxyz(2, il) = rxyz(2, i) - alat(2)
4674 rxyz(3, il) = rxyz(3, i) + alat(3)
4675 END DO
4676
4677 in = icell(0, l1, ll2 - 1, ll3 - 1)
4678 icell(0, l1, -1, -1) = in
4679 DO ii = 1, in
4680 i = icell(ii, l1, ll2 - 1, ll3 - 1)
4681 il = il + 1
4682 IF (il > nn) cpabort("enlarge laymx")
4683 lay(il) = i
4684 icell(ii, l1, -1, -1) = il
4685 rxyz(1, il) = rxyz(1, i)
4686 rxyz(2, il) = rxyz(2, i) - alat(2)
4687 rxyz(3, il) = rxyz(3, i) - alat(3)
4688 END DO
4689
4690 END DO
4691
4692! y axis
4693 DO l2 = 0, ll2 - 1
4694
4695 in = icell(0, 0, l2, 0)
4696 icell(0, ll1, l2, ll3) = in
4697 DO ii = 1, in
4698 i = icell(ii, 0, l2, 0)
4699 il = il + 1
4700 IF (il > nn) cpabort("enlarge laymx")
4701 lay(il) = i
4702 icell(ii, ll1, l2, ll3) = il
4703 rxyz(1, il) = rxyz(1, i) + alat(1)
4704 rxyz(2, il) = rxyz(2, i)
4705 rxyz(3, il) = rxyz(3, i) + alat(3)
4706 END DO
4707
4708 in = icell(0, 0, l2, ll3 - 1)
4709 icell(0, ll1, l2, -1) = in
4710 DO ii = 1, in
4711 i = icell(ii, 0, l2, ll3 - 1)
4712 il = il + 1
4713 IF (il > nn) cpabort("enlarge laymx")
4714 lay(il) = i
4715 icell(ii, ll1, l2, -1) = il
4716 rxyz(1, il) = rxyz(1, i) + alat(1)
4717 rxyz(2, il) = rxyz(2, i)
4718 rxyz(3, il) = rxyz(3, i) - alat(3)
4719 END DO
4720
4721 in = icell(0, ll1 - 1, l2, 0)
4722 icell(0, -1, l2, ll3) = in
4723 DO ii = 1, in
4724 i = icell(ii, ll1 - 1, l2, 0)
4725 il = il + 1
4726 IF (il > nn) cpabort("enlarge laymx")
4727 lay(il) = i
4728 icell(ii, -1, l2, ll3) = il
4729 rxyz(1, il) = rxyz(1, i) - alat(1)
4730 rxyz(2, il) = rxyz(2, i)
4731 rxyz(3, il) = rxyz(3, i) + alat(3)
4732 END DO
4733
4734 in = icell(0, ll1 - 1, l2, ll3 - 1)
4735 icell(0, -1, l2, -1) = in
4736 DO ii = 1, in
4737 i = icell(ii, ll1 - 1, l2, ll3 - 1)
4738 il = il + 1
4739 IF (il > nn) cpabort("enlarge laymx")
4740 lay(il) = i
4741 icell(ii, -1, l2, -1) = il
4742 rxyz(1, il) = rxyz(1, i) - alat(1)
4743 rxyz(2, il) = rxyz(2, i)
4744 rxyz(3, il) = rxyz(3, i) - alat(3)
4745 END DO
4746
4747 END DO
4748
4749! z axis
4750 DO l3 = 0, ll3 - 1
4751
4752 in = icell(0, 0, 0, l3)
4753 icell(0, ll1, ll2, l3) = in
4754 DO ii = 1, in
4755 i = icell(ii, 0, 0, l3)
4756 il = il + 1
4757 IF (il > nn) cpabort("enlarge laymx")
4758 lay(il) = i
4759 icell(ii, ll1, ll2, l3) = il
4760 rxyz(1, il) = rxyz(1, i) + alat(1)
4761 rxyz(2, il) = rxyz(2, i) + alat(2)
4762 rxyz(3, il) = rxyz(3, i)
4763 END DO
4764
4765 in = icell(0, ll1 - 1, 0, l3)
4766 icell(0, -1, ll2, l3) = in
4767 DO ii = 1, in
4768 i = icell(ii, ll1 - 1, 0, l3)
4769 il = il + 1
4770 IF (il > nn) cpabort("enlarge laymx")
4771 lay(il) = i
4772 icell(ii, -1, ll2, l3) = il
4773 rxyz(1, il) = rxyz(1, i) - alat(1)
4774 rxyz(2, il) = rxyz(2, i) + alat(2)
4775 rxyz(3, il) = rxyz(3, i)
4776 END DO
4777
4778 in = icell(0, 0, ll2 - 1, l3)
4779 icell(0, ll1, -1, l3) = in
4780 DO ii = 1, in
4781 i = icell(ii, 0, ll2 - 1, l3)
4782 il = il + 1
4783 IF (il > nn) cpabort("enlarge laymx")
4784 lay(il) = i
4785 icell(ii, ll1, -1, l3) = il
4786 rxyz(1, il) = rxyz(1, i) + alat(1)
4787 rxyz(2, il) = rxyz(2, i) - alat(2)
4788 rxyz(3, il) = rxyz(3, i)
4789 END DO
4790
4791 in = icell(0, ll1 - 1, ll2 - 1, l3)
4792 icell(0, -1, -1, l3) = in
4793 DO ii = 1, in
4794 i = icell(ii, ll1 - 1, ll2 - 1, l3)
4795 il = il + 1
4796 IF (il > nn) cpabort("enlarge laymx")
4797 lay(il) = i
4798 icell(ii, -1, -1, l3) = il
4799 rxyz(1, il) = rxyz(1, i) - alat(1)
4800 rxyz(2, il) = rxyz(2, i) - alat(2)
4801 rxyz(3, il) = rxyz(3, i)
4802 END DO
4803
4804 END DO
4805
4806! corners
4807 in = icell(0, 0, 0, 0)
4808 icell(0, ll1, ll2, ll3) = in
4809 DO ii = 1, in
4810 i = icell(ii, 0, 0, 0)
4811 il = il + 1
4812 IF (il > nn) cpabort("enlarge laymx")
4813 lay(il) = i
4814 icell(ii, ll1, ll2, ll3) = il
4815 rxyz(1, il) = rxyz(1, i) + alat(1)
4816 rxyz(2, il) = rxyz(2, i) + alat(2)
4817 rxyz(3, il) = rxyz(3, i) + alat(3)
4818 END DO
4819
4820 in = icell(0, ll1 - 1, 0, 0)
4821 icell(0, -1, ll2, ll3) = in
4822 DO ii = 1, in
4823 i = icell(ii, ll1 - 1, 0, 0)
4824 il = il + 1
4825 IF (il > nn) cpabort("enlarge laymx")
4826 lay(il) = i
4827 icell(ii, -1, ll2, ll3) = il
4828 rxyz(1, il) = rxyz(1, i) - alat(1)
4829 rxyz(2, il) = rxyz(2, i) + alat(2)
4830 rxyz(3, il) = rxyz(3, i) + alat(3)
4831 END DO
4832
4833 in = icell(0, 0, ll2 - 1, 0)
4834 icell(0, ll1, -1, ll3) = in
4835 DO ii = 1, in
4836 i = icell(ii, 0, ll2 - 1, 0)
4837 il = il + 1
4838 IF (il > nn) cpabort("enlarge laymx")
4839 lay(il) = i
4840 icell(ii, ll1, -1, ll3) = il
4841 rxyz(1, il) = rxyz(1, i) + alat(1)
4842 rxyz(2, il) = rxyz(2, i) - alat(2)
4843 rxyz(3, il) = rxyz(3, i) + alat(3)
4844 END DO
4845
4846 in = icell(0, ll1 - 1, ll2 - 1, 0)
4847 icell(0, -1, -1, ll3) = in
4848 DO ii = 1, in
4849 i = icell(ii, ll1 - 1, ll2 - 1, 0)
4850 il = il + 1
4851 IF (il > nn) cpabort("enlarge laymx")
4852 lay(il) = i
4853 icell(ii, -1, -1, ll3) = il
4854 rxyz(1, il) = rxyz(1, i) - alat(1)
4855 rxyz(2, il) = rxyz(2, i) - alat(2)
4856 rxyz(3, il) = rxyz(3, i) + alat(3)
4857 END DO
4858
4859 in = icell(0, 0, 0, ll3 - 1)
4860 icell(0, ll1, ll2, -1) = in
4861 DO ii = 1, in
4862 i = icell(ii, 0, 0, ll3 - 1)
4863 il = il + 1
4864 IF (il > nn) cpabort("enlarge laymx")
4865 lay(il) = i
4866 icell(ii, ll1, ll2, -1) = il
4867 rxyz(1, il) = rxyz(1, i) + alat(1)
4868 rxyz(2, il) = rxyz(2, i) + alat(2)
4869 rxyz(3, il) = rxyz(3, i) - alat(3)
4870 END DO
4871
4872 in = icell(0, ll1 - 1, 0, ll3 - 1)
4873 icell(0, -1, ll2, -1) = in
4874 DO ii = 1, in
4875 i = icell(ii, ll1 - 1, 0, ll3 - 1)
4876 il = il + 1
4877 IF (il > nn) cpabort("enlarge laymx")
4878 lay(il) = i
4879 icell(ii, -1, ll2, -1) = il
4880 rxyz(1, il) = rxyz(1, i) - alat(1)
4881 rxyz(2, il) = rxyz(2, i) + alat(2)
4882 rxyz(3, il) = rxyz(3, i) - alat(3)
4883 END DO
4884
4885 in = icell(0, 0, ll2 - 1, ll3 - 1)
4886 icell(0, ll1, -1, -1) = in
4887 DO ii = 1, in
4888 i = icell(ii, 0, ll2 - 1, ll3 - 1)
4889 il = il + 1
4890 IF (il > nn) cpabort("enlarge laymx")
4891 lay(il) = i
4892 icell(ii, ll1, -1, -1) = il
4893 rxyz(1, il) = rxyz(1, i) + alat(1)
4894 rxyz(2, il) = rxyz(2, i) - alat(2)
4895 rxyz(3, il) = rxyz(3, i) - alat(3)
4896 END DO
4897
4898 in = icell(0, ll1 - 1, ll2 - 1, ll3 - 1)
4899 icell(0, -1, -1, -1) = in
4900 DO ii = 1, in
4901 i = icell(ii, ll1 - 1, ll2 - 1, ll3 - 1)
4902 il = il + 1
4903 IF (il > nn) cpabort("enlarge laymx")
4904 lay(il) = i
4905 icell(ii, -1, -1, -1) = il
4906 rxyz(1, il) = rxyz(1, i) - alat(1)
4907 rxyz(2, il) = rxyz(2, i) - alat(2)
4908 rxyz(3, il) = rxyz(3, i) - alat(3)
4909 END DO
4910
4911 nnbrx = 3*nnbrx/2
4912 ALLOCATE (rel(5, nnbrx*nat))
4913 indlstx = 0
4914
4915 npr = 1
4916 iam = 0
4917
4918 cut2 = cut**2
4919! assign contiguous portions of the arrays lstb and rel to the threads
4920 myspace = (nat*nnbrx)/npr
4921 IF (iam == 0) myspaceout = myspace
4922! Verlet list, relative positions
4923 indlst = 0
4924 DO l3 = 0, ll3 - 1
4925 DO l2 = 0, ll2 - 1
4926 DO l1 = 0, ll1 - 1
4927 DO ii = 1, icell(0, l1, l2, l3)
4928 iat = icell(ii, l1, l2, l3)
4929 IF (((iat - 1)*npr)/nat == iam) THEN
4930! write(6,*) 'sublstiat:iam,iat',iam,iat
4931 lsta(1, iat) = iam*myspace + indlst + 1
4932 CALL tersoff_sublstiat_l(iat, nn, ncx, ll1, ll2, ll3, l1, l2, l3, myspace, &
4933 rxyz, icell, lstb(iam*myspace + 1), lay, rel(1, iam*myspace + 1), cut2, indlst)
4934 lsta(2, iat) = iam*myspace + indlst
4935 ipb = lsta(1, iat)
4936 ndat = lsta(2, iat) - lsta(1, iat) + 1
4937 END IF
4938
4939 END DO
4940 END DO
4941 END DO
4942 END DO
4943 indlstx = max(indlstx, indlst)
4944
4945 IF (indlstx >= myspaceout) cpabort("NNBRX too small")
4946 npr = 1
4947 iam = 0
4948
4949 npjx = 300; npjkx = 6000
4950 istopg = 0
4951!end of creating pairlist part------------------------------------------------------------
4952!Energy-----------------------------------------------------------------------------------
4953 urtot = 0.0d0
4954 nptot = 0
4955
4956 f = 0.0d0
4957 uatot = 0.0d0
4958 do_i: DO i = 1, nmol
4959 CALL tersoff_subeniat_l(i, nmol, npmax, kinds, x, r1, r2, cr, ca, alr, ala, xyzrrefdf, uadurdf, urtot, lsta, lstb, nnbrx, rel, pi)
4960 END DO do_i
4961
4962 urtot = 0.5d0*urtot
4963!Force------------------------------------------------------------------------------------
4964 f = 0.0d0
4965 uatot = 0.0d0
4966
4967 do_if: DO i = 1, nmol
4968 CALL tersoff_subfiat_l(i,nmol,npmax,nnmax,kinds,pn,co_bcd,bcsq,dsq,h,xyzrrefdf,uadurdf,f,uatot,dkeij,lsta,lstb,nnbrx)
4969 END DO do_if
4970
4971 f = 0.5d0*f
4972 uatot = 0.5d0*uatot
4973!-----------------------------------------------------------------------------------------
4974 DEALLOCATE (rxyz, icell, lay, rel)
4975 RETURN
4976 END SUBROUTINE tersoff_pairlist_energy_forces
4977
4978!-----------------------------------------------------------------------------------------
4979! **************************************************************************************************
4980!> \brief ...
4981!> \param iat ...
4982!> \param nn ...
4983!> \param ncx ...
4984!> \param ll1 ...
4985!> \param ll2 ...
4986!> \param ll3 ...
4987!> \param l1 ...
4988!> \param l2 ...
4989!> \param l3 ...
4990!> \param myspace ...
4991!> \param rxyz ...
4992!> \param icell ...
4993!> \param lstb ...
4994!> \param lay ...
4995!> \param rel ...
4996!> \param cut2 ...
4997!> \param indlst ...
4998! **************************************************************************************************
4999 SUBROUTINE tersoff_sublstiat_l(iat, nn, ncx, ll1, ll2, ll3, l1, l2, l3, myspace, &
5000 rxyz, icell, lstb, lay, rel, cut2, indlst)
5001! finds the neighbours of atom iat (specified by lsta and lstb) and and
5002! the relative position rel of iat with respect to these neighbours
5003 INTEGER :: iat, nn, ncx, ll1, ll2, ll3, l1, l2, l3, &
5004 myspace
5005 REAL(8) :: rxyz(3, nn)
5006 INTEGER :: icell(0:ncx, -1:ll1, -1:ll2, -1:ll3), lstb(0:myspace - 1), lay(nn)
5007 REAL(8) :: rel(5, 0:myspace - 1), cut2
5008 INTEGER :: indlst
5009
5010 INTEGER :: jat, jj, k1, k2, k3
5011 REAL(8) :: rr2, tt, tti, xrel, yrel, zrel
5012
5013 DO k3 = l3 - 1, l3 + 1
5014 DO k2 = l2 - 1, l2 + 1
5015 DO k1 = l1 - 1, l1 + 1
5016 DO jj = 1, icell(0, k1, k2, k3)
5017 jat = icell(jj, k1, k2, k3)
5018 IF (jat == iat) cycle
5019 xrel = rxyz(1, iat) - rxyz(1, jat)
5020 yrel = rxyz(2, iat) - rxyz(2, jat)
5021 zrel = rxyz(3, iat) - rxyz(3, jat)
5022 rr2 = xrel**2 + yrel**2 + zrel**2
5023 IF (rr2 <= cut2) THEN
5024 indlst = min(indlst, myspace - 1)
5025 lstb(indlst) = lay(jat)
5026! write(6,*) 'iat,indlst,lay(jat)',iat,indlst,lay(jat)
5027 tt = sqrt(rr2)
5028 tti = 1.d0/tt
5029 rel(1, indlst) = xrel*tti
5030 rel(2, indlst) = yrel*tti
5031 rel(3, indlst) = zrel*tti
5032 rel(4, indlst) = tt
5033 rel(5, indlst) = tti
5034 indlst = indlst + 1
5035 END IF
5036 END DO
5037 END DO
5038 END DO
5039 END DO
5040
5041 RETURN
5042 END SUBROUTINE tersoff_sublstiat_l
5043
5044! **************************************************************************************************
5045!> \brief ...
5046!> \param i ...
5047!> \param Nmol ...
5048!> \param Npmax ...
5049!> \param Kinds ...
5050!> \param X ...
5051!> \param R1 ...
5052!> \param R2 ...
5053!> \param Cr ...
5054!> \param Ca ...
5055!> \param alr ...
5056!> \param ala ...
5057!> \param XYZRrefdf ...
5058!> \param UadUrdf ...
5059!> \param Urtot ...
5060!> \param lsta ...
5061!> \param lstb ...
5062!> \param nnbrx ...
5063!> \param rel ...
5064!> \param pi ...
5065! **************************************************************************************************
5066 SUBROUTINE tersoff_subeniat_l(i,Nmol,Npmax,Kinds,X,R1,R2,Cr,Ca,alr,ala,XYZRrefdf,UadUrdf,Urtot,lsta,lstb,nnbrx,rel,pi)
5067 INTEGER :: i
5068 INTEGER, INTENT(in) :: nmol, npmax
5069 INTEGER, DIMENSION(1:Nmol), INTENT(in) :: kinds
5070 REAL(8), DIMENSION(1:2, 1:2), INTENT(in) :: x, r1, r2, cr, ca, alr, ala
5071 REAL(8), DIMENSION(1:6*Npmax), INTENT(inout) :: xyzrrefdf
5072 REAL(8), DIMENSION(1:3*Npmax), INTENT(inout) :: uadurdf
5073 REAL(8), INTENT(inout) :: urtot
5074 INTEGER, INTENT(in) :: lsta(2, nmol), nnbrx, lstb(nnbrx*nmol)
5075 REAL(8), INTENT(in) :: rel(5, nnbrx*nmol)
5076 REAL(8) :: pi
5077
5078 INTEGER :: j, ki, kj, l, nppt3, nppt6, nptot
5079 REAL(8) :: alaij, alrij, dfij, fij, pl1, pl2, r1ij, &
5080 r2ij, rij, rreij, ua, ur, xij, yij, zij
5081
5082! #######################################
5083! # Calculate XYZRrefdf, UadUrdf, Urtot #
5084! #######################################
5085 ki = kinds(i)
5086
5087 do_j: DO l = lsta(1, i), lsta(2, i)
5088 j = lstb(l)
5089
5090 kj = kinds(j)
5091 r2ij = r2(ki, kj)
5092 rij = rel(4, l)
5093 xij = rel(1, l)
5094 yij = rel(2, l)
5095 zij = rel(3, l)
5096 nptot = l
5097
5098 nppt3 = 3*(nptot - 1)
5099 nppt6 = 6*(nptot - 1)
5100 rreij = rel(5, l)
5101
5102 xyzrrefdf(nppt6 + 1) = xij
5103 xyzrrefdf(nppt6 + 2) = yij
5104 xyzrrefdf(nppt6 + 3) = zij
5105 xyzrrefdf(nppt6 + 4) = rreij
5106
5107 alrij = alr(ki, kj)
5108 alaij = ala(ki, kj)
5109
5110 ur = cr(ki, kj)*dexp(-alrij*rij)
5111 ua = -ca(ki, kj)*dexp(-alaij*rij)*x(ki, kj)
5112 r1ij = r1(ki, kj)
5113
5114 IF (rij <= r1ij) THEN
5115 xyzrrefdf(nppt6 + 5) = 1.0d0
5116 xyzrrefdf(nppt6 + 6) = 0.0d0
5117 urtot = urtot + ur
5118 uadurdf(nppt3 + 1) = ua
5119 uadurdf(nppt3 + 2) = -alrij*ur
5120 uadurdf(nppt3 + 3) = -alaij*ua
5121 ELSE
5122 pl1 = pi/(r2ij - r1ij)
5123 pl2 = pl1*(rij - r1ij)
5124 fij = 0.5d0 + 0.5d0*dcos(pl2)
5125 dfij = -0.5d0*pl1*dsin(pl2)
5126 xyzrrefdf(nppt6 + 5) = fij
5127 xyzrrefdf(nppt6 + 6) = dfij
5128 urtot = urtot + fij*ur
5129 uadurdf(nppt3 + 1) = fij*ua
5130 uadurdf(nppt3 + 2) = (dfij - alrij*fij)*ur
5131 uadurdf(nppt3 + 3) = (dfij - alaij*fij)*ua
5132 END IF
5133 END DO do_j
5134 END SUBROUTINE tersoff_subeniat_l
5135
5136! **************************************************************************************************
5137!> \brief ...
5138!> \param i ...
5139!> \param Nmol ...
5140!> \param Npmax ...
5141!> \param NNmax ...
5142!> \param Kinds ...
5143!> \param Pn ...
5144!> \param Co_bcd ...
5145!> \param bcsq ...
5146!> \param dsq ...
5147!> \param h ...
5148!> \param XYZRrefdf ...
5149!> \param UadUrdf ...
5150!> \param F ...
5151!> \param Uatot ...
5152!> \param dkEij ...
5153!> \param lsta ...
5154!> \param lstb ...
5155!> \param nnbrx ...
5156! **************************************************************************************************
5157 SUBROUTINE tersoff_subfiat_l(i,Nmol,Npmax,NNmax,Kinds,Pn,Co_bcd,bcsq,dsq,h,XYZRrefdf,UadUrdf,F,Uatot,dkEij,lsta,lstb,nnbrx)
5158 INTEGER :: i
5159 INTEGER, INTENT(in) :: nmol, npmax, nnmax
5160 INTEGER, DIMENSION(1:Nmol), INTENT(in) :: kinds
5161 REAL(8), DIMENSION(1:2), INTENT(in) :: pn, co_bcd, bcsq, dsq, h
5162 REAL(8), DIMENSION(1:6*Npmax), INTENT(in) :: xyzrrefdf
5163 REAL(8), DIMENSION(1:3*Npmax), INTENT(in) :: uadurdf
5164 REAL(8), DIMENSION(1:3*Nmol), INTENT(inout) :: f
5165 REAL(8), INTENT(inout) :: uatot
5166 REAL(8), DIMENSION(1:3*NNmax) :: dkeij
5167 INTEGER, INTENT(in) :: lsta(2, nmol), nnbrx, lstb(nnbrx*nmol)
5168
5169 INTEGER :: ij, ijpt3, ijpt6, ik, ikpt6, ipb, ipe, &
5170 ipt3, jpt3, ki, kpt3, nkpt3
5171 REAL(8) :: bcsqi, bij, co1_dkeij, co2_dkeij, co_cdi, co_dhcosi, co_hcosi, co_mb1, co_mb2, &
5172 co_pa, cosijk, dfij, dfik, dfxi, dfxj, dfxk, dfyi, dfyj, dfyk, dfzi, dfzj, dfzk, dgi, &
5173 djeij, dsqi, dxjeij2, dyjeij2, dzjeij2, eij, fdg, fdgcos, fij, fik, gi, hi, pni, rreij, &
5174 rreik, ua, xrreij, xrreik, yrreij, yrreik, zrreij, zrreik
5175
5176 ipb = lsta(1, i)
5177 ipe = lsta(2, i)
5178
5179 ki = kinds(i)
5180 bcsqi = bcsq(ki)
5181 dsqi = dsq(ki)
5182 hi = h(ki)
5183 pni = pn(ki)
5184
5185 co_cdi = co_bcd(ki)
5186
5187 dfxi = 0.0d0
5188 dfyi = 0.0d0
5189 dfzi = 0.0d0
5190
5191 do_j: DO ij = ipb, ipe, +1
5192
5193 ijpt3 = 3*(ij - 1)
5194 ijpt6 = 6*(ij - 1)
5195
5196 xrreij = xyzrrefdf(ijpt6 + 1)
5197 yrreij = xyzrrefdf(ijpt6 + 2)
5198 zrreij = xyzrrefdf(ijpt6 + 3)
5199 rreij = xyzrrefdf(ijpt6 + 4)
5200 fij = xyzrrefdf(ijpt6 + 5)
5201 dfij = xyzrrefdf(ijpt6 + 6)
5202
5203 eij = 0.0d0
5204 djeij = 0.0d0
5205 dxjeij2 = 0.0d0
5206 dyjeij2 = 0.0d0
5207 dzjeij2 = 0.0d0
5208
5209 nkpt3 = -3
5210 do_k: DO ik = ipb, ipe, +1
5211
5212 nkpt3 = nkpt3 + 3
5213
5214 ikij: IF (ik /= ij) THEN
5215
5216 ikpt6 = 6*(ik - 1)
5217
5218 xrreik = xyzrrefdf(ikpt6 + 1)
5219 yrreik = xyzrrefdf(ikpt6 + 2)
5220 zrreik = xyzrrefdf(ikpt6 + 3)
5221 rreik = xyzrrefdf(ikpt6 + 4)
5222 fik = xyzrrefdf(ikpt6 + 5)
5223 dfik = xyzrrefdf(ikpt6 + 6)
5224
5225 cosijk = xrreij*xrreik + yrreij*yrreik + zrreij*zrreik
5226
5227 co_hcosi = hi - cosijk
5228 co_dhcosi = 1.0d0/(dsqi + co_hcosi*co_hcosi)
5229 gi = -bcsqi*co_dhcosi
5230 dgi = 2.0d0*co_hcosi*co_dhcosi*gi
5231 gi = gi + co_cdi
5232
5233 eij = eij + fik*gi
5234
5235 fdg = fik*dgi
5236 fdgcos = fdg*cosijk
5237
5238 djeij = djeij + fdgcos
5239
5240 dxjeij2 = dxjeij2 + fdg*xrreik
5241 dyjeij2 = dyjeij2 + fdg*yrreik
5242 dzjeij2 = dzjeij2 + fdg*zrreik
5243
5244 co1_dkeij = -dfik*gi + fdgcos*rreik
5245 co2_dkeij = -fdg*rreik
5246
5247 dkeij(nkpt3 + 1) = co1_dkeij*xrreik + co2_dkeij*xrreij
5248 dkeij(nkpt3 + 2) = co1_dkeij*yrreik + co2_dkeij*yrreij
5249 dkeij(nkpt3 + 3) = co1_dkeij*zrreik + co2_dkeij*zrreij
5250
5251 ELSE
5252 dkeij(nkpt3 + 1) = 0.0d0
5253 dkeij(nkpt3 + 2) = 0.0d0
5254 dkeij(nkpt3 + 3) = 0.0d0
5255 END IF ikij
5256
5257 END DO do_k
5258
5259 bij = 1.0d0 + eij**pni
5260 ua = uadurdf(ijpt3 + 1)*bij**(-0.5d0/pni)
5261 uatot = uatot + ua
5262
5263 co_pa = uadurdf(ijpt3 + 2) + uadurdf(ijpt3 + 3)*bij**(-0.5d0/pni)
5264
5265 ceij: IF (nkpt3 > 0) THEN
5266
5267 co_mb1 = ua*0.5d0*eij**(pni - 1.0d0)/bij
5268 co_mb2 = co_mb1*rreij
5269
5270 nkpt3 = -3
5271 DO ik = ipb, ipe, +1
5272
5273 nkpt3 = nkpt3 + 3
5274 dfxk = co_mb1*dkeij(nkpt3 + 1)
5275 dfyk = co_mb1*dkeij(nkpt3 + 2)
5276 dfzk = co_mb1*dkeij(nkpt3 + 3)
5277
5278 kpt3 = 3*(lstb(ik) - 1)
5279 f(kpt3 + 1) = f(kpt3 + 1) + dfxk
5280 f(kpt3 + 2) = f(kpt3 + 2) + dfyk
5281 f(kpt3 + 3) = f(kpt3 + 3) + dfzk
5282
5283 dfxi = dfxi + dfxk
5284 dfyi = dfyi + dfyk
5285 dfzi = dfzi + dfzk
5286
5287 END DO
5288
5289 dfxj = co_pa*xrreij + co_mb2*(xrreij*djeij - dxjeij2)
5290 dfyj = co_pa*yrreij + co_mb2*(yrreij*djeij - dyjeij2)
5291 dfzj = co_pa*zrreij + co_mb2*(zrreij*djeij - dzjeij2)
5292
5293 ELSE
5294
5295 dfxj = co_pa*xrreij
5296 dfyj = co_pa*yrreij
5297 dfzj = co_pa*zrreij
5298
5299 END IF ceij
5300
5301 jpt3 = 3*(lstb(ij) - 1)
5302
5303 f(jpt3 + 1) = f(jpt3 + 1) + dfxj
5304 f(jpt3 + 2) = f(jpt3 + 2) + dfyj
5305 f(jpt3 + 3) = f(jpt3 + 3) + dfzj
5306
5307 dfxi = dfxi + dfxj
5308 dfyi = dfyi + dfyj
5309 dfzi = dfzi + dfzj
5310
5311 END DO do_j
5312
5313 ipt3 = 3*(i - 1)
5314 f(ipt3 + 1) = f(ipt3 + 1) - dfxi
5315 f(ipt3 + 2) = f(ipt3 + 2) - dfyi
5316 f(ipt3 + 3) = f(ipt3 + 3) - dfzi
5317 END SUBROUTINE tersoff_subfiat_l
5318
5319END MODULE eip_silicon
modulo
static GRID_HOST_DEVICE int modulo(int a, int m)
Equivalent of Fortran's MODULO, which always return a positive number. https://gcc....
Definition grid_common.h:125
atomic_kind_list_types
represent a simple array based list of the given type
Definition atomic_kind_list_types.F:15
atomic_kind_types
Define the atomic kind types and their sub types.
Definition atomic_kind_types.F:23
atomic_kind_types::get_atomic_kind
subroutine, public get_atomic_kind(atomic_kind, fist_potential, element_symbol, name, mass, kind_number, natom, atom_list, rcov, rvdw, z, qeff, apol, cpol, mm_radius, shell, shell_active, damping)
Get attributes of an atomic kind.
Definition atomic_kind_types.F:138
cell_types
Handles all functions related to the CELL.
Definition cell_types.F:15
cell_types::get_cell
subroutine, public get_cell(cell, alpha, beta, gamma, deth, orthorhombic, abc, periodic, h, h_inv, symmetry_id, tag)
Get informations about a simulation cell.
Definition cell_types.F:200
cp_log_handling
various routines to log and control the output. The idea is that decisions about where to log should ...
Definition cp_log_handling.F:41
cp_log_handling::cp_get_default_logger
type(cp_logger_type) function, pointer, public cp_get_default_logger()
returns the default logger
Definition cp_log_handling.F:234
cp_output_handling
routines to handle the output, The idea is to remove the decision of wheter to output and what to out...
Definition cp_output_handling.F:25
cp_output_handling::cp_print_key_unit_nr
integer function, public cp_print_key_unit_nr(logger, basis_section, print_key_path, extension, middle_name, local, log_filename, ignore_should_output, file_form, file_position, file_action, file_status, do_backup, on_file, is_new_file, mpi_io, fout)
...
Definition cp_output_handling.F:853
cp_output_handling::cp_print_key_finished_output
subroutine, public cp_print_key_finished_output(unit_nr, logger, basis_section, print_key_path, local, ignore_should_output, on_file, mpi_io)
should be called after you finish working with a unit obtained with cp_print_key_unit_nr,...
Definition cp_output_handling.F:1063
cp_output_handling::cp_p_file
integer, parameter, public cp_p_file
Definition cp_output_handling.F:103
cp_output_handling::cp_print_key_should_output
integer function, public cp_print_key_should_output(iteration_info, basis_section, print_key_path, used_print_key, first_time)
returns what should be done with the given property if btest(res,cp_p_store) then the property should...
Definition cp_output_handling.F:345
cp_subsys_types
types that represent a subsys, i.e. a part of the system
Definition cp_subsys_types.F:16
cp_subsys_types::cp_subsys_get
subroutine, public cp_subsys_get(subsys, ref_count, atomic_kinds, atomic_kind_set, particles, particle_set, local_particles, molecules, molecule_set, molecule_kinds, molecule_kind_set, local_molecules, para_env, colvar_p, shell_particles, core_particles, gci, multipoles, natom, nparticle, ncore, nshell, nkind, atprop, virial, results, cell, cell_ref, use_ref_cell)
returns information about various attributes of the given subsys
Definition cp_subsys_types.F:364
distribution_1d_types
stores a lists of integer that are local to a processor. The idea is that these integers represent ob...
Definition distribution_1d_types.F:24
eip_environment_types
The environment for the empirical interatomic potential methods.
Definition eip_environment_types.F:14
eip_environment_types::eip_env_get
subroutine, public eip_env_get(eip_env, eip_model, eip_energy, eip_energy_var, eip_forces, coord_avg, coord_var, count, subsys, atomic_kind_set, particle_set, local_particles, molecule_kind_set, molecule_set, local_molecules, eip_input, force_env_input, cell, cell_ref, use_ref_cell, eip_kinetic_energy, eip_potential_energy, virial)
Returns various attributes of the eip environment.
Definition eip_environment_types.F:174
eip_silicon
Empirical interatomic potentials for Silicon.
Definition eip_silicon.F:17
eip_silicon::eip_stillinger_weber
subroutine, public eip_stillinger_weber(eip_env)
Interface routine of the Stillinger-Weber force field to CP2K.
Definition eip_silicon.F:405
eip_silicon::eip_lenosky
subroutine, public eip_lenosky(eip_env)
Interface routine of Goedecker's Lenosky force field to CP2K.
Definition eip_silicon.F:241
eip_silicon::eip_tersoff
subroutine, public eip_tersoff(eip_env)
Interface routine of the Tersoff force field to CP2K.
Definition eip_silicon.F:567
eip_silicon::eip_bazant
subroutine, public eip_bazant(eip_env)
Interface routine of Goedecker's Bazant EDIP to CP2K.
Definition eip_silicon.F:75
input_section_types
objects that represent the structure of input sections and the data contained in an input section
Definition input_section_types.F:15
input_section_types::section_vals_get_subs_vals
recursive type(section_vals_type) function, pointer, public section_vals_get_subs_vals(section_vals, subsection_name, i_rep_section, can_return_null)
returns the values of the requested subsection
Definition input_section_types.F:731
kinds
Defines the basic variable types.
Definition kinds.F:23
kinds::dp
integer, parameter, public dp
Definition kinds.F:34
message_passing
Interface to the message passing library MPI.
Definition message_passing.F:23
particle_types
Define the data structure for the particle information.
Definition particle_types.F:19
physcon
Definition of physical constants:
Definition physcon.F:68
physcon::evolt
real(kind=dp), parameter, public evolt
Definition physcon.F:183
physcon::angstrom
real(kind=dp), parameter, public angstrom
Definition physcon.F:144
atomic_kind_list_types::atomic_kind_list_type
represent a list of objects
Definition atomic_kind_list_types.F:45
atomic_kind_types::atomic_kind_type
Provides all information about an atomic kind.
Definition atomic_kind_types.F:45
cell_types::cell_type
Type defining parameters related to the simulation cell.
Definition cell_types.F:60
cp_log_handling::cp_logger_type
type of a logger, at the moment it contains just a print level starting at which level it should be l...
Definition cp_log_handling.F:140
cp_subsys_types::cp_subsys_type
represents a system: atoms, molecules, their pos,vel,...
Definition cp_subsys_types.F:89
distribution_1d_types::distribution_1d_type
structure to store local (to a processor) ordered lists of integers.
Definition distribution_1d_types.F:62
eip_environment_types::eip_environment_type
The empirical interatomic potential environment.
Definition eip_environment_types.F:83
input_section_types::section_vals_type
stores the values of a section
Definition input_section_types.F:127
message_passing::mp_para_env_type
stores all the informations relevant to an mpi environment
Definition message_passing.F:743
particle_types::particle_type
Definition particle_types.F:35