(git:b279b6b)
eip_silicon.F
Go to the documentation of this file.
1 !--------------------------------------------------------------------------------------------------!
2 ! CP2K: A general program to perform molecular dynamics simulations !
3 ! Copyright 2000-2024 CP2K developers group <https://cp2k.org> !
4 ! !
5 ! SPDX-License-Identifier: GPL-2.0-or-later !
6 !--------------------------------------------------------------------------------------------------!
7 
8 ! **************************************************************************************************
9 !> \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@phys.chem.ethz.ch)
16 ! **************************************************************************************************
17 MODULE 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
52 
53 !***
54 
55 CONTAINS
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@phys.chem.ethz.ch)
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@phys.chem.ethz.ch)
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 Print routine for the EIP energies
396 !> \param eip_env The eip environment of matter
397 !> \param output_unit The output unit
398 !> \par History
399 !> 03.2006 initial create [tdk]
400 !> \author Thomas D. Kuehne (tkuehne@phys.chem.ethz.ch)
401 !> \note
402 !> As usual the EIP energies differ from the DFT energies!
403 !> Only the relative energy differences are correctly reproduced.
404 ! **************************************************************************************************
405  SUBROUTINE eip_print_energies(eip_env, output_unit)
406  TYPE(eip_environment_type), POINTER :: eip_env
407  INTEGER, INTENT(IN) :: output_unit
408 
409 ! ------------------------------------------------------------------------
410 
411  IF (output_unit > 0) THEN
412  WRITE (unit=output_unit, fmt="(/,(T3,A,T55,F25.14))") &
413  "Kinetic energy [Hartree]: ", eip_env%eip_kinetic_energy, &
414  "Potential energy [Hartree]: ", eip_env%eip_potential_energy, &
415  "Total EIP energy [Hartree]: ", eip_env%eip_energy
416  END IF
417 
418  END SUBROUTINE eip_print_energies
419 
420 ! **************************************************************************************************
421 !> \brief Print routine for the variance of the energy/atom
422 !> \param eip_env The eip environment of matter
423 !> \param output_unit The output unit
424 !> \par History
425 !> 03.2006 initial create [tdk]
426 !> \author Thomas D. Kuehne (tkuehne@phys.chem.ethz.ch)
427 ! **************************************************************************************************
428  SUBROUTINE eip_print_energy_var(eip_env, output_unit)
429  TYPE(eip_environment_type), POINTER :: eip_env
430  INTEGER, INTENT(IN) :: output_unit
431 
432  INTEGER :: unit_nr
433 
434 ! ------------------------------------------------------------------------
435 
436  unit_nr = output_unit
437 
438  IF (unit_nr > 0) THEN
439 
440  WRITE (unit_nr, *) ""
441  WRITE (unit_nr, *) "The variance of the EIP energy/atom!"
442  WRITE (unit_nr, *) ""
443  WRITE (unit_nr, *) eip_env%eip_energy_var
444 
445  END IF
446 
447  END SUBROUTINE eip_print_energy_var
448 
449 ! **************************************************************************************************
450 !> \brief Print routine for the forces
451 !> \param eip_env The eip environment of matter
452 !> \param output_unit The output unit
453 !> \par History
454 !> 03.2006 initial create [tdk]
455 !> \author Thomas D. Kuehne (tkuehne@phys.chem.ethz.ch)
456 ! **************************************************************************************************
457  SUBROUTINE eip_print_forces(eip_env, output_unit)
458  TYPE(eip_environment_type), POINTER :: eip_env
459  INTEGER, INTENT(IN) :: output_unit
460 
461  INTEGER :: iatom, natom, unit_nr
462  TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
463 
464 ! ------------------------------------------------------------------------
465 
466  NULLIFY (particle_set)
467 
468  unit_nr = output_unit
469 
470  IF (unit_nr > 0) THEN
471 
472  CALL eip_env_get(eip_env=eip_env, particle_set=particle_set)
473 
474  natom = SIZE(particle_set)
475 
476  WRITE (unit_nr, *) ""
477  WRITE (unit_nr, *) "The EIP forces!"
478  WRITE (unit_nr, *) ""
479  WRITE (unit_nr, *) "Total EIP forces [Hartree/Bohr]"
480  DO iatom = 1, natom
481  WRITE (unit_nr, *) eip_env%eip_forces(1:3, iatom)
482  END DO
483 
484  END IF
485 
486  END SUBROUTINE eip_print_forces
487 
488 ! **************************************************************************************************
489 !> \brief Print routine for the average coordination number
490 !> \param eip_env The eip environment of matter
491 !> \param output_unit The output unit
492 !> \par History
493 !> 03.2006 initial create [tdk]
494 !> \author Thomas D. Kuehne (tkuehne@phys.chem.ethz.ch)
495 ! **************************************************************************************************
496  SUBROUTINE eip_print_coord_avg(eip_env, output_unit)
497  TYPE(eip_environment_type), POINTER :: eip_env
498  INTEGER, INTENT(IN) :: output_unit
499 
500  INTEGER :: unit_nr
501 
502 ! ------------------------------------------------------------------------
503 
504  unit_nr = output_unit
505 
506  IF (unit_nr > 0) THEN
507 
508  WRITE (unit_nr, *) ""
509  WRITE (unit_nr, *) "The average coordination number!"
510  WRITE (unit_nr, *) ""
511  WRITE (unit_nr, *) eip_env%coord_avg
512 
513  END IF
514 
515  END SUBROUTINE eip_print_coord_avg
516 
517 ! **************************************************************************************************
518 !> \brief Print routine for the variance of the coordination number
519 !> \param eip_env The eip environment of matter
520 !> \param output_unit The output unit
521 !> \par History
522 !> 03.2006 initial create [tdk]
523 !> \author Thomas D. Kuehne (tkuehne@phys.chem.ethz.ch)
524 ! **************************************************************************************************
525  SUBROUTINE eip_print_coord_var(eip_env, output_unit)
526  TYPE(eip_environment_type), POINTER :: eip_env
527  INTEGER, INTENT(IN) :: output_unit
528 
529  INTEGER :: unit_nr
530 
531 ! ------------------------------------------------------------------------
532 
533  unit_nr = output_unit
534 
535  IF (unit_nr > 0) THEN
536 
537  WRITE (unit_nr, *) ""
538  WRITE (unit_nr, *) "The variance of the coordination number!"
539  WRITE (unit_nr, *) ""
540  WRITE (unit_nr, *) eip_env%coord_var
541 
542  END IF
543 
544  END SUBROUTINE eip_print_coord_var
545 
546 ! **************************************************************************************************
547 !> \brief Print routine for the function call counter
548 !> \param eip_env The eip environment of matter
549 !> \param output_unit The output unit
550 !> \par History
551 !> 03.2006 initial create [tdk]
552 !> \author Thomas D. Kuehne (tkuehne@phys.chem.ethz.ch)
553 ! **************************************************************************************************
554  SUBROUTINE eip_print_count(eip_env, output_unit)
555  TYPE(eip_environment_type), POINTER :: eip_env
556  INTEGER, INTENT(IN) :: output_unit
557 
558  INTEGER :: unit_nr
559 
560 ! ------------------------------------------------------------------------
561 
562  unit_nr = output_unit
563 
564  IF (unit_nr > 0) THEN
565 
566  WRITE (unit_nr, *) ""
567  WRITE (unit_nr, *) "The function call counter!"
568  WRITE (unit_nr, *) ""
569  WRITE (unit_nr, *) eip_env%count
570 
571  END IF
572 
573  END SUBROUTINE eip_print_count
574 
575 ! **************************************************************************************************
576 !> \brief Bazant's EDIP (environment-dependent interatomic potential) for Silicon
577 !> by Stefan Goedecker
578 !> \param nat number of atoms
579 !> \param alat lattice constants of the orthorombic box containing the particles
580 !> \param rxyz0 atomic positions in Angstrom, may be modified on output.
581 !> If an atom is outside the box the program will bring it back
582 !> into the box by translations through alat
583 !> \param fxyz forces in eV/A
584 !> \param ener total energy in eV
585 !> \param coord average coordination number
586 !> \param ener_var variance of the energy/atom
587 !> \param coord_var variance of the coordination number
588 !> \param count count is increased by one per call, has to be initialized
589 !> to 0.e0_dp before first call of eip_bazant
590 !> \par Literature
591 !> http://www-math.mit.edu/~bazant/EDIP
592 !> M.Z. Bazant & E. Kaxiras: Modeling of Covalent Bonding in Solids by
593 !> Inversion of Cohesive Energy Curves;
594 !> Phys. Rev. Lett. 77, 4370 (1996)
595 !> M.Z. Bazant, E. Kaxiras and J.F. Justo: Environment-dependent interatomic
596 !> potential for bulk silicon;
597 !> Phys. Rev. B 56, 8542-8552 (1997)
598 !> S. Goedecker: Optimization and parallelization of a force field for silicon
599 !> using OpenMP; CPC 148, 1 (2002)
600 !> \par History
601 !> 03.2006 initial create [tdk]
602 !> \author Thomas D. Kuehne (tkuehne@phys.chem.ethz.ch)
603 ! **************************************************************************************************
604  SUBROUTINE eip_bazant_silicon(nat, alat, rxyz0, fxyz, ener, coord, ener_var, &
605  coord_var, count)
606 
607  INTEGER :: nat
608  REAL(kind=dp) :: alat, rxyz0, fxyz, ener, coord, &
609  ener_var, coord_var, count
610 
611  dimension rxyz0(3, nat), fxyz(3, nat), alat(3)
612  REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :) :: rxyz
613  INTEGER, ALLOCATABLE, DIMENSION(:, :) :: lsta
614  INTEGER, ALLOCATABLE, DIMENSION(:) :: lstb
615  INTEGER, ALLOCATABLE, DIMENSION(:) :: lay
616  INTEGER, ALLOCATABLE, DIMENSION(:, :, :, :) :: icell
617  REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :) :: rel
618  REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :) :: txyz
619  REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :) :: s2, s3, sz
620  INTEGER, ALLOCATABLE, DIMENSION(:) :: num2, num3, numz
621 
622  REAL(kind=dp) :: coord2, cut, cut2, ener2, rlc1i, rlc2i, rlc3i, tcoord, &
623  tcoord2, tener, tener2
624  INTEGER :: iam, iat, iat1, iat2, ii, i, il, in, indlst, indlstx, istop, &
625  istopg, l2, l3, laymx, ll1, ll2, ll3, lot, max_nbrs, myspace, &
626  l1, myspaceout, ncx, nn, nnbrx, npr
627 
628 ! cut=par_a
629  cut = 3.1213820e0_dp + 1.e-14_dp
630 
631  IF (count .EQ. 0) OPEN (unit=10, file='bazant.mon', status='unknown')
632  count = count + 1.e0_dp
633 
634 ! linear scaling calculation of verlet list
635  ll1 = int(alat(1)/cut)
636  IF (ll1 .LT. 1) cpabort("alat(1) too small")
637  ll2 = int(alat(2)/cut)
638  IF (ll2 .LT. 1) cpabort("alat(2) too small")
639  ll3 = int(alat(3)/cut)
640  IF (ll3 .LT. 1) cpabort("alat(3) too small")
641 
642 ! determine number of threads
643  npr = 1
644 !$OMP PARALLEL PRIVATE(iam) SHARED (npr) DEFAULT(NONE)
645 !$ iam = omp_get_thread_num()
646 !$ if (iam .eq. 0) npr = omp_get_num_threads()
647 !$OMP END PARALLEL
648 
649 ! linear scaling calculation of verlet list
650 
651  IF (npr .LE. 1) THEN !serial if too few processors to gain by parallelizing
652 
653 ! set ncx for serial case, ncx for parallel case set below
654  ncx = 16
655  loop_ncx_s: DO
656  ALLOCATE (icell(0:ncx, -1:ll1, -1:ll2, -1:ll3))
657  icell(0, -1:ll1, -1:ll2, -1:ll3) = 0
658  rlc1i = ll1/alat(1)
659  rlc2i = ll2/alat(2)
660  rlc3i = ll3/alat(3)
661 
662  loop_iat_s: DO iat = 1, nat
663  rxyz0(1, iat) = modulo(modulo(rxyz0(1, iat), alat(1)), alat(1))
664  rxyz0(2, iat) = modulo(modulo(rxyz0(2, iat), alat(2)), alat(2))
665  rxyz0(3, iat) = modulo(modulo(rxyz0(3, iat), alat(3)), alat(3))
666  l1 = int(rxyz0(1, iat)*rlc1i)
667  l2 = int(rxyz0(2, iat)*rlc2i)
668  l3 = int(rxyz0(3, iat)*rlc3i)
669 
670  ii = icell(0, l1, l2, l3)
671  ii = ii + 1
672  icell(0, l1, l2, l3) = ii
673  IF (ii .GT. ncx) THEN
674  WRITE (10, *) count, 'NCX too small', ncx
675  DEALLOCATE (icell)
676  ncx = ncx*2
677  cycle loop_ncx_s
678  END IF
679  icell(ii, l1, l2, l3) = iat
680  END DO loop_iat_s
681  EXIT loop_ncx_s
682  END DO loop_ncx_s
683 
684  ELSE ! parallel case
685 
686 ! periodization of particles can be done in parallel
687 !$OMP PARALLEL DO SHARED (alat,nat,rxyz0) PRIVATE(iat) DEFAULT(NONE)
688  DO iat = 1, nat
689  rxyz0(1, iat) = modulo(modulo(rxyz0(1, iat), alat(1)), alat(1))
690  rxyz0(2, iat) = modulo(modulo(rxyz0(2, iat), alat(2)), alat(2))
691  rxyz0(3, iat) = modulo(modulo(rxyz0(3, iat), alat(3)), alat(3))
692  END DO
693 !$OMP END PARALLEL DO
694 
695 ! assignment to cell is done serially
696 ! set ncx for parallel case, ncx for serial case set above
697  ncx = 16
698  loop_ncx_p: DO
699  ALLOCATE (icell(0:ncx, -1:ll1, -1:ll2, -1:ll3))
700  icell(0, -1:ll1, -1:ll2, -1:ll3) = 0
701 
702  rlc1i = ll1/alat(1)
703  rlc2i = ll2/alat(2)
704  rlc3i = ll3/alat(3)
705 
706  loop_iat_p: DO iat = 1, nat
707  l1 = int(rxyz0(1, iat)*rlc1i)
708  l2 = int(rxyz0(2, iat)*rlc2i)
709  l3 = int(rxyz0(3, iat)*rlc3i)
710  ii = icell(0, l1, l2, l3)
711  ii = ii + 1
712  icell(0, l1, l2, l3) = ii
713  IF (ii .GT. ncx) THEN
714  WRITE (10, *) count, 'NCX too small', ncx
715  DEALLOCATE (icell)
716  ncx = ncx*2
717  cycle loop_ncx_p
718  END IF
719  icell(ii, l1, l2, l3) = iat
720  END DO loop_iat_p
721  EXIT loop_ncx_p
722  END DO loop_ncx_p
723 
724  END IF
725 
726 ! duplicate all atoms within boundary layer
727  laymx = ncx*(2*ll1*ll2 + 2*ll1*ll3 + 2*ll2*ll3 + 4*ll1 + 4*ll2 + 4*ll3 + 8)
728  nn = nat + laymx
729  ALLOCATE (rxyz(3, nn), lay(nn))
730  DO iat = 1, nat
731  lay(iat) = iat
732  rxyz(1, iat) = rxyz0(1, iat)
733  rxyz(2, iat) = rxyz0(2, iat)
734  rxyz(3, iat) = rxyz0(3, iat)
735  END DO
736  il = nat
737 ! xy plane
738  DO l2 = 0, ll2 - 1
739  DO l1 = 0, ll1 - 1
740 
741  in = icell(0, l1, l2, 0)
742  icell(0, l1, l2, ll3) = in
743  DO ii = 1, in
744  i = icell(ii, l1, l2, 0)
745  il = il + 1
746  IF (il .GT. nn) cpabort("enlarge laymx")
747  lay(il) = i
748  icell(ii, l1, l2, ll3) = il
749  rxyz(1, il) = rxyz(1, i)
750  rxyz(2, il) = rxyz(2, i)
751  rxyz(3, il) = rxyz(3, i) + alat(3)
752  END DO
753 
754  in = icell(0, l1, l2, ll3 - 1)
755  icell(0, l1, l2, -1) = in
756  DO ii = 1, in
757  i = icell(ii, l1, l2, ll3 - 1)
758  il = il + 1
759  IF (il .GT. nn) cpabort("enlarge laymx")
760  lay(il) = i
761  icell(ii, l1, l2, -1) = il
762  rxyz(1, il) = rxyz(1, i)
763  rxyz(2, il) = rxyz(2, i)
764  rxyz(3, il) = rxyz(3, i) - alat(3)
765  END DO
766 
767  END DO
768  END DO
769 
770 ! yz plane
771  DO l3 = 0, ll3 - 1
772  DO l2 = 0, ll2 - 1
773 
774  in = icell(0, 0, l2, l3)
775  icell(0, ll1, l2, l3) = in
776  DO ii = 1, in
777  i = icell(ii, 0, l2, l3)
778  il = il + 1
779  IF (il .GT. nn) cpabort("enlarge laymx")
780  lay(il) = i
781  icell(ii, ll1, l2, l3) = il
782  rxyz(1, il) = rxyz(1, i) + alat(1)
783  rxyz(2, il) = rxyz(2, i)
784  rxyz(3, il) = rxyz(3, i)
785  END DO
786 
787  in = icell(0, ll1 - 1, l2, l3)
788  icell(0, -1, l2, l3) = in
789  DO ii = 1, in
790  i = icell(ii, ll1 - 1, l2, l3)
791  il = il + 1
792  IF (il .GT. nn) cpabort("enlarge laymx")
793  lay(il) = i
794  icell(ii, -1, l2, l3) = il
795  rxyz(1, il) = rxyz(1, i) - alat(1)
796  rxyz(2, il) = rxyz(2, i)
797  rxyz(3, il) = rxyz(3, i)
798  END DO
799 
800  END DO
801  END DO
802 
803 ! xz plane
804  DO l3 = 0, ll3 - 1
805  DO l1 = 0, ll1 - 1
806 
807  in = icell(0, l1, 0, l3)
808  icell(0, l1, ll2, l3) = in
809  DO ii = 1, in
810  i = icell(ii, l1, 0, l3)
811  il = il + 1
812  IF (il .GT. nn) cpabort("enlarge laymx")
813  lay(il) = i
814  icell(ii, l1, ll2, l3) = il
815  rxyz(1, il) = rxyz(1, i)
816  rxyz(2, il) = rxyz(2, i) + alat(2)
817  rxyz(3, il) = rxyz(3, i)
818  END DO
819 
820  in = icell(0, l1, ll2 - 1, l3)
821  icell(0, l1, -1, l3) = in
822  DO ii = 1, in
823  i = icell(ii, l1, ll2 - 1, l3)
824  il = il + 1
825  IF (il .GT. nn) cpabort("enlarge laymx")
826  lay(il) = i
827  icell(ii, l1, -1, l3) = il
828  rxyz(1, il) = rxyz(1, i)
829  rxyz(2, il) = rxyz(2, i) - alat(2)
830  rxyz(3, il) = rxyz(3, i)
831  END DO
832 
833  END DO
834  END DO
835 
836 ! x axis
837  DO l1 = 0, ll1 - 1
838 
839  in = icell(0, l1, 0, 0)
840  icell(0, l1, ll2, ll3) = in
841  DO ii = 1, in
842  i = icell(ii, l1, 0, 0)
843  il = il + 1
844  IF (il .GT. nn) cpabort("enlarge laymx")
845  lay(il) = i
846  icell(ii, l1, ll2, ll3) = il
847  rxyz(1, il) = rxyz(1, i)
848  rxyz(2, il) = rxyz(2, i) + alat(2)
849  rxyz(3, il) = rxyz(3, i) + alat(3)
850  END DO
851 
852  in = icell(0, l1, 0, ll3 - 1)
853  icell(0, l1, ll2, -1) = in
854  DO ii = 1, in
855  i = icell(ii, l1, 0, ll3 - 1)
856  il = il + 1
857  IF (il .GT. nn) cpabort("enlarge laymx")
858  lay(il) = i
859  icell(ii, l1, ll2, -1) = il
860  rxyz(1, il) = rxyz(1, i)
861  rxyz(2, il) = rxyz(2, i) + alat(2)
862  rxyz(3, il) = rxyz(3, i) - alat(3)
863  END DO
864 
865  in = icell(0, l1, ll2 - 1, 0)
866  icell(0, l1, -1, ll3) = in
867  DO ii = 1, in
868  i = icell(ii, l1, ll2 - 1, 0)
869  il = il + 1
870  IF (il .GT. nn) cpabort("enlarge laymx")
871  lay(il) = i
872  icell(ii, l1, -1, ll3) = il
873  rxyz(1, il) = rxyz(1, i)
874  rxyz(2, il) = rxyz(2, i) - alat(2)
875  rxyz(3, il) = rxyz(3, i) + alat(3)
876  END DO
877 
878  in = icell(0, l1, ll2 - 1, ll3 - 1)
879  icell(0, l1, -1, -1) = in
880  DO ii = 1, in
881  i = icell(ii, l1, ll2 - 1, ll3 - 1)
882  il = il + 1
883  IF (il .GT. nn) cpabort("enlarge laymx")
884  lay(il) = i
885  icell(ii, l1, -1, -1) = il
886  rxyz(1, il) = rxyz(1, i)
887  rxyz(2, il) = rxyz(2, i) - alat(2)
888  rxyz(3, il) = rxyz(3, i) - alat(3)
889  END DO
890 
891  END DO
892 
893 ! y axis
894  DO l2 = 0, ll2 - 1
895 
896  in = icell(0, 0, l2, 0)
897  icell(0, ll1, l2, ll3) = in
898  DO ii = 1, in
899  i = icell(ii, 0, l2, 0)
900  il = il + 1
901  IF (il .GT. nn) cpabort("enlarge laymx")
902  lay(il) = i
903  icell(ii, ll1, l2, ll3) = il
904  rxyz(1, il) = rxyz(1, i) + alat(1)
905  rxyz(2, il) = rxyz(2, i)
906  rxyz(3, il) = rxyz(3, i) + alat(3)
907  END DO
908 
909  in = icell(0, 0, l2, ll3 - 1)
910  icell(0, ll1, l2, -1) = in
911  DO ii = 1, in
912  i = icell(ii, 0, l2, ll3 - 1)
913  il = il + 1
914  IF (il .GT. nn) cpabort("enlarge laymx")
915  lay(il) = i
916  icell(ii, ll1, l2, -1) = il
917  rxyz(1, il) = rxyz(1, i) + alat(1)
918  rxyz(2, il) = rxyz(2, i)
919  rxyz(3, il) = rxyz(3, i) - alat(3)
920  END DO
921 
922  in = icell(0, ll1 - 1, l2, 0)
923  icell(0, -1, l2, ll3) = in
924  DO ii = 1, in
925  i = icell(ii, ll1 - 1, l2, 0)
926  il = il + 1
927  IF (il .GT. nn) cpabort("enlarge laymx")
928  lay(il) = i
929  icell(ii, -1, l2, ll3) = il
930  rxyz(1, il) = rxyz(1, i) - alat(1)
931  rxyz(2, il) = rxyz(2, i)
932  rxyz(3, il) = rxyz(3, i) + alat(3)
933  END DO
934 
935  in = icell(0, ll1 - 1, l2, ll3 - 1)
936  icell(0, -1, l2, -1) = in
937  DO ii = 1, in
938  i = icell(ii, ll1 - 1, l2, ll3 - 1)
939  il = il + 1
940  IF (il .GT. nn) cpabort("enlarge laymx")
941  lay(il) = i
942  icell(ii, -1, l2, -1) = il
943  rxyz(1, il) = rxyz(1, i) - alat(1)
944  rxyz(2, il) = rxyz(2, i)
945  rxyz(3, il) = rxyz(3, i) - alat(3)
946  END DO
947 
948  END DO
949 
950 ! z axis
951  DO l3 = 0, ll3 - 1
952 
953  in = icell(0, 0, 0, l3)
954  icell(0, ll1, ll2, l3) = in
955  DO ii = 1, in
956  i = icell(ii, 0, 0, l3)
957  il = il + 1
958  IF (il .GT. nn) cpabort("enlarge laymx")
959  lay(il) = i
960  icell(ii, ll1, ll2, l3) = il
961  rxyz(1, il) = rxyz(1, i) + alat(1)
962  rxyz(2, il) = rxyz(2, i) + alat(2)
963  rxyz(3, il) = rxyz(3, i)
964  END DO
965 
966  in = icell(0, ll1 - 1, 0, l3)
967  icell(0, -1, ll2, l3) = in
968  DO ii = 1, in
969  i = icell(ii, ll1 - 1, 0, l3)
970  il = il + 1
971  IF (il .GT. nn) cpabort("enlarge laymx")
972  lay(il) = i
973  icell(ii, -1, ll2, l3) = il
974  rxyz(1, il) = rxyz(1, i) - alat(1)
975  rxyz(2, il) = rxyz(2, i) + alat(2)
976  rxyz(3, il) = rxyz(3, i)
977  END DO
978 
979  in = icell(0, 0, ll2 - 1, l3)
980  icell(0, ll1, -1, l3) = in
981  DO ii = 1, in
982  i = icell(ii, 0, ll2 - 1, l3)
983  il = il + 1
984  IF (il .GT. nn) cpabort("enlarge laymx")
985  lay(il) = i
986  icell(ii, ll1, -1, l3) = il
987  rxyz(1, il) = rxyz(1, i) + alat(1)
988  rxyz(2, il) = rxyz(2, i) - alat(2)
989  rxyz(3, il) = rxyz(3, i)
990  END DO
991 
992  in = icell(0, ll1 - 1, ll2 - 1, l3)
993  icell(0, -1, -1, l3) = in
994  DO ii = 1, in
995  i = icell(ii, ll1 - 1, ll2 - 1, l3)
996  il = il + 1
997  IF (il .GT. nn) cpabort("enlarge laymx")
998  lay(il) = i
999  icell(ii, -1, -1, l3) = il
1000  rxyz(1, il) = rxyz(1, i) - alat(1)
1001  rxyz(2, il) = rxyz(2, i) - alat(2)
1002  rxyz(3, il) = rxyz(3, i)
1003  END DO
1004 
1005  END DO
1006 
1007 ! corners
1008  in = icell(0, 0, 0, 0)
1009  icell(0, ll1, ll2, ll3) = in
1010  DO ii = 1, in
1011  i = icell(ii, 0, 0, 0)
1012  il = il + 1
1013  IF (il .GT. nn) cpabort("enlarge laymx")
1014  lay(il) = i
1015  icell(ii, ll1, ll2, ll3) = il
1016  rxyz(1, il) = rxyz(1, i) + alat(1)
1017  rxyz(2, il) = rxyz(2, i) + alat(2)
1018  rxyz(3, il) = rxyz(3, i) + alat(3)
1019  END DO
1020 
1021  in = icell(0, ll1 - 1, 0, 0)
1022  icell(0, -1, ll2, ll3) = in
1023  DO ii = 1, in
1024  i = icell(ii, ll1 - 1, 0, 0)
1025  il = il + 1
1026  IF (il .GT. nn) cpabort("enlarge laymx")
1027  lay(il) = i
1028  icell(ii, -1, ll2, ll3) = il
1029  rxyz(1, il) = rxyz(1, i) - alat(1)
1030  rxyz(2, il) = rxyz(2, i) + alat(2)
1031  rxyz(3, il) = rxyz(3, i) + alat(3)
1032  END DO
1033 
1034  in = icell(0, 0, ll2 - 1, 0)
1035  icell(0, ll1, -1, ll3) = in
1036  DO ii = 1, in
1037  i = icell(ii, 0, ll2 - 1, 0)
1038  il = il + 1
1039  IF (il .GT. nn) cpabort("enlarge laymx")
1040  lay(il) = i
1041  icell(ii, ll1, -1, ll3) = il
1042  rxyz(1, il) = rxyz(1, i) + alat(1)
1043  rxyz(2, il) = rxyz(2, i) - alat(2)
1044  rxyz(3, il) = rxyz(3, i) + alat(3)
1045  END DO
1046 
1047  in = icell(0, ll1 - 1, ll2 - 1, 0)
1048  icell(0, -1, -1, ll3) = in
1049  DO ii = 1, in
1050  i = icell(ii, ll1 - 1, ll2 - 1, 0)
1051  il = il + 1
1052  IF (il .GT. nn) cpabort("enlarge laymx")
1053  lay(il) = i
1054  icell(ii, -1, -1, ll3) = il
1055  rxyz(1, il) = rxyz(1, i) - alat(1)
1056  rxyz(2, il) = rxyz(2, i) - alat(2)
1057  rxyz(3, il) = rxyz(3, i) + alat(3)
1058  END DO
1059 
1060  in = icell(0, 0, 0, ll3 - 1)
1061  icell(0, ll1, ll2, -1) = in
1062  DO ii = 1, in
1063  i = icell(ii, 0, 0, ll3 - 1)
1064  il = il + 1
1065  IF (il .GT. nn) cpabort("enlarge laymx")
1066  lay(il) = i
1067  icell(ii, ll1, ll2, -1) = il
1068  rxyz(1, il) = rxyz(1, i) + alat(1)
1069  rxyz(2, il) = rxyz(2, i) + alat(2)
1070  rxyz(3, il) = rxyz(3, i) - alat(3)
1071  END DO
1072 
1073  in = icell(0, ll1 - 1, 0, ll3 - 1)
1074  icell(0, -1, ll2, -1) = in
1075  DO ii = 1, in
1076  i = icell(ii, ll1 - 1, 0, ll3 - 1)
1077  il = il + 1
1078  IF (il .GT. nn) cpabort("enlarge laymx")
1079  lay(il) = i
1080  icell(ii, -1, ll2, -1) = il
1081  rxyz(1, il) = rxyz(1, i) - alat(1)
1082  rxyz(2, il) = rxyz(2, i) + alat(2)
1083  rxyz(3, il) = rxyz(3, i) - alat(3)
1084  END DO
1085 
1086  in = icell(0, 0, ll2 - 1, ll3 - 1)
1087  icell(0, ll1, -1, -1) = in
1088  DO ii = 1, in
1089  i = icell(ii, 0, ll2 - 1, ll3 - 1)
1090  il = il + 1
1091  IF (il .GT. nn) cpabort("enlarge laymx")
1092  lay(il) = i
1093  icell(ii, ll1, -1, -1) = il
1094  rxyz(1, il) = rxyz(1, i) + alat(1)
1095  rxyz(2, il) = rxyz(2, i) - alat(2)
1096  rxyz(3, il) = rxyz(3, i) - alat(3)
1097  END DO
1098 
1099  in = icell(0, ll1 - 1, ll2 - 1, ll3 - 1)
1100  icell(0, -1, -1, -1) = in
1101  DO ii = 1, in
1102  i = icell(ii, ll1 - 1, ll2 - 1, ll3 - 1)
1103  il = il + 1
1104  IF (il .GT. nn) cpabort("enlarge laymx")
1105  lay(il) = i
1106  icell(ii, -1, -1, -1) = il
1107  rxyz(1, il) = rxyz(1, i) - alat(1)
1108  rxyz(2, il) = rxyz(2, i) - alat(2)
1109  rxyz(3, il) = rxyz(3, i) - alat(3)
1110  END DO
1111 
1112  ALLOCATE (lsta(2, nat))
1113  nnbrx = 12
1114  loop_nnbrx: DO
1115  ALLOCATE (lstb(nnbrx*nat), rel(5, nnbrx*nat))
1116 
1117  indlstx = 0
1118 
1119 !$OMP PARALLEL DEFAULT(NONE) &
1120 !$OMP PRIVATE(iat,cut2,iam,ii,indlst,l1,l2,l3,myspace,npr) &
1121 !$OMP SHARED (indlstx,nat,nn,nnbrx,ncx,ll1,ll2,ll3,icell,lsta,lstb,lay, &
1122 !$OMP rel,rxyz,cut,myspaceout)
1123 
1124  npr = 1
1125 !$ npr = omp_get_num_threads()
1126  iam = 0
1127 !$ iam = omp_get_thread_num()
1128 
1129  cut2 = cut**2
1130 ! assign contiguous portions of the arrays lstb and rel to the threads
1131  myspace = (nat*nnbrx)/npr
1132  IF (iam .EQ. 0) myspaceout = myspace
1133 ! Verlet list, relative positions
1134  indlst = 0
1135  loop_l3: DO l3 = 0, ll3 - 1
1136  loop_l2: DO l2 = 0, ll2 - 1
1137  loop_l1: DO l1 = 0, ll1 - 1
1138  loop_ii: DO ii = 1, icell(0, l1, l2, l3)
1139  iat = icell(ii, l1, l2, l3)
1140  IF (((iat - 1)*npr)/nat .EQ. iam) THEN
1141 ! write(*,*) 'sublstiat:iam,iat',iam,iat
1142  lsta(1, iat) = iam*myspace + indlst + 1
1143  CALL sublstiat_b(iat, nn, ncx, ll1, ll2, ll3, l1, l2, l3, myspace, &
1144  rxyz, icell, lstb(iam*myspace + 1), lay, &
1145  rel(1, iam*myspace + 1), cut2, indlst)
1146  lsta(2, iat) = iam*myspace + indlst
1147 ! write(*,'(a,4(x,i3),100(x,i2))') &
1148 ! 'iam,iat,lsta',iam,iat,lsta(1,iat),lsta(2,iat), &
1149 ! (lstb(j),j=lsta(1,iat),lsta(2,iat))
1150  END IF
1151  END DO loop_ii
1152  END DO loop_l1
1153  END DO loop_l2
1154  END DO loop_l3
1155 !$OMP CRITICAL
1156  indlstx = max(indlstx, indlst)
1157 !$OMP END CRITICAL
1158 !$OMP END PARALLEL
1159 
1160  IF (indlstx .GE. myspaceout) THEN
1161  WRITE (10, *) count, 'NNBRX too small', nnbrx
1162  DEALLOCATE (lstb, rel)
1163  nnbrx = 3*nnbrx/2
1164  cycle loop_nnbrx
1165  END IF
1166  EXIT loop_nnbrx
1167  END DO loop_nnbrx
1168 
1169  istopg = 0
1170 
1171 !$OMP PARALLEL DEFAULT(NONE) &
1172 !$OMP PRIVATE(iam,npr,iat,iat1,iat2,lot,istop,tcoord,tcoord2, &
1173 !$OMP tener,tener2,txyz,s2,s3,sz,num2,num3,numz,max_nbrs) &
1174 !$OMP SHARED (nat,nnbrx,lsta,lstb,rel,ener,ener2,fxyz,coord,coord2,istopg)
1175 
1176  npr = 1
1177 !$ npr = omp_get_num_threads()
1178  iam = 0
1179 !$ iam = omp_get_thread_num()
1180 
1181  max_nbrs = 30
1182 
1183  IF (npr .NE. 1) THEN
1184 ! PARALLEL CASE
1185 ! create temporary private scalars for reduction sum on energies and
1186 ! temporary private array for reduction sum on forces
1187 !$OMP CRITICAL
1188  ALLOCATE (txyz(3, nat), s2(max_nbrs, 8), s3(max_nbrs, 7), sz(max_nbrs, 6), &
1189  num2(max_nbrs), num3(max_nbrs), numz(max_nbrs))
1190 !$OMP END CRITICAL
1191  IF (iam .EQ. 0) THEN
1192  ener = 0.e0_dp
1193  ener2 = 0.e0_dp
1194  coord = 0.e0_dp
1195  coord2 = 0.e0_dp
1196  END IF
1197 !$OMP DO
1198  DO iat = 1, nat
1199  fxyz(1, iat) = 0.e0_dp
1200  fxyz(2, iat) = 0.e0_dp
1201  fxyz(3, iat) = 0.e0_dp
1202  END DO
1203 !$OMP BARRIER
1204 
1205 ! Each thread treats at most lot atoms
1206  lot = int(real(nat, kind=dp)/real(npr, kind=dp) + .999999999999e0_dp)
1207  iat1 = iam*lot + 1
1208  iat2 = min((iam + 1)*lot, nat)
1209 ! write(*,*) 'subfeniat:iat1,iat2,iam',iat1,iat2,iam
1210  CALL subfeniat_b(iat1, iat2, nat, lsta, lstb, rel, tener, tener2, &
1211  tcoord, tcoord2, nnbrx, txyz, max_nbrs, istop, &
1212  s2(1, 1), s2(1, 2), s2(1, 3), s2(1, 4), s2(1, 5), s2(1, 6), s2(1, 7), s2(1, 8), &
1213  num2, s3(1, 1), s3(1, 2), s3(1, 3), s3(1, 4), s3(1, 5), s3(1, 6), s3(1, 7), &
1214  num3, sz(1, 1), sz(1, 2), sz(1, 3), sz(1, 4), sz(1, 5), sz(1, 6), numz)
1215 
1216 !$OMP CRITICAL
1217  ener = ener + tener
1218  ener2 = ener2 + tener2
1219  coord = coord + tcoord
1220  coord2 = coord2 + tcoord2
1221  istopg = istopg + istop
1222  DO iat = 1, nat
1223  fxyz(1, iat) = fxyz(1, iat) + txyz(1, iat)
1224  fxyz(2, iat) = fxyz(2, iat) + txyz(2, iat)
1225  fxyz(3, iat) = fxyz(3, iat) + txyz(3, iat)
1226  END DO
1227  DEALLOCATE (txyz, s2, s3, sz, num2, num3, numz)
1228 !$OMP END CRITICAL
1229 
1230  ELSE
1231 ! SERIAL CASE
1232  iat1 = 1
1233  iat2 = nat
1234  ALLOCATE (s2(max_nbrs, 8), s3(max_nbrs, 7), sz(max_nbrs, 6), &
1235  num2(max_nbrs), num3(max_nbrs), numz(max_nbrs))
1236  CALL subfeniat_b(iat1, iat2, nat, lsta, lstb, rel, ener, ener2, &
1237  coord, coord2, nnbrx, fxyz, max_nbrs, istopg, &
1238  s2(1, 1), s2(1, 2), s2(1, 3), s2(1, 4), s2(1, 5), s2(1, 6), s2(1, 7), s2(1, 8), &
1239  num2, s3(1, 1), s3(1, 2), s3(1, 3), s3(1, 4), s3(1, 5), s3(1, 6), s3(1, 7), &
1240  num3, sz(1, 1), sz(1, 2), sz(1, 3), sz(1, 4), sz(1, 5), sz(1, 6), numz)
1241  DEALLOCATE (s2, s3, sz, num2, num3, numz)
1242 
1243  END IF
1244 !$OMP END PARALLEL
1245 
1246 ! write(*,*) 'ener,norm force', &
1247 ! ener,DNRM2(3*nat,fxyz,1)
1248  IF (istopg .GT. 0) cpabort("DIMENSION ERROR (see WARNING above)")
1249  ener_var = ener2/nat - (ener/nat)**2
1250  coord = coord/nat
1251  coord_var = coord2/nat - coord**2
1252 
1253  DEALLOCATE (rxyz, icell, lay, lsta, lstb, rel)
1254 
1255  END SUBROUTINE eip_bazant_silicon
1256 
1257 ! **************************************************************************************************
1258 !> \brief ...
1259 !> \param iat1 ...
1260 !> \param iat2 ...
1261 !> \param nat ...
1262 !> \param lsta ...
1263 !> \param lstb ...
1264 !> \param rel ...
1265 !> \param ener ...
1266 !> \param ener2 ...
1267 !> \param coord ...
1268 !> \param coord2 ...
1269 !> \param nnbrx ...
1270 !> \param ff ...
1271 !> \param max_nbrs ...
1272 !> \param istop ...
1273 !> \param s2_t0 ...
1274 !> \param s2_t1 ...
1275 !> \param s2_t2 ...
1276 !> \param s2_t3 ...
1277 !> \param s2_dx ...
1278 !> \param s2_dy ...
1279 !> \param s2_dz ...
1280 !> \param s2_r ...
1281 !> \param num2 ...
1282 !> \param s3_g ...
1283 !> \param s3_dg ...
1284 !> \param s3_rinv ...
1285 !> \param s3_dx ...
1286 !> \param s3_dy ...
1287 !> \param s3_dz ...
1288 !> \param s3_r ...
1289 !> \param num3 ...
1290 !> \param sz_df ...
1291 !> \param sz_sum ...
1292 !> \param sz_dx ...
1293 !> \param sz_dy ...
1294 !> \param sz_dz ...
1295 !> \param sz_r ...
1296 !> \param numz ...
1297 ! **************************************************************************************************
1298  SUBROUTINE subfeniat_b(iat1, iat2, nat, lsta, lstb, rel, ener, ener2, &
1299  coord, coord2, nnbrx, ff, max_nbrs, istop, &
1300  s2_t0, s2_t1, s2_t2, s2_t3, s2_dx, s2_dy, s2_dz, s2_r, &
1301  num2, s3_g, s3_dg, s3_rinv, s3_dx, s3_dy, s3_dz, s3_r, &
1302  num3, sz_df, sz_sum, sz_dx, sz_dy, sz_dz, sz_r, numz)
1303 ! This subroutine is a modification of a subroutine that is available at
1304 ! http://www-math.mit.edu/~bazant/EDIP/ and for which Martin Z. Bazant
1305 ! and Harvard University have a 1997 copyright.
1306 ! The modifications were done by S. Goedecker on April 10, 2002.
1307 ! The routines are included with the permission of M. Bazant into this package.
1308 
1309 ! ------------------------- VARIABLE DECLARATIONS -------------------------
1310  INTEGER :: iat1, iat2, nat, lsta(2, nat)
1311  REAL(kind=dp) :: ener, ener2, coord, coord2
1312  INTEGER :: nnbrx
1313  REAL(kind=dp) :: rel(5, nnbrx*nat)
1314  INTEGER :: lstb(nnbrx*nat)
1315  REAL(kind=dp) :: ff(3, nat)
1316  INTEGER :: max_nbrs, istop
1317  REAL(kind=dp) :: s2_t0(max_nbrs), s2_t1(max_nbrs), s2_t2(max_nbrs), s2_t3(max_nbrs), &
1318  s2_dx(max_nbrs), s2_dy(max_nbrs), s2_dz(max_nbrs), s2_r(max_nbrs)
1319  INTEGER :: num2(max_nbrs)
1320  REAL(kind=dp) :: s3_g(max_nbrs), s3_dg(max_nbrs), s3_rinv(max_nbrs), s3_dx(max_nbrs), &
1321  s3_dy(max_nbrs), s3_dz(max_nbrs), s3_r(max_nbrs)
1322  INTEGER :: num3(max_nbrs)
1323  REAL(kind=dp) :: sz_df(max_nbrs), sz_sum(max_nbrs), &
1324  sz_dx(max_nbrs), sz_dy(max_nbrs), &
1325  sz_dz(max_nbrs), sz_r(max_nbrs)
1326  INTEGER :: numz(max_nbrs)
1327 
1328  INTEGER :: i, j, k, l, n, n2, n3, nj, nk, nl, nz
1329  REAL(kind=dp) :: bmc, cmbinv, coord_iat, dedrl, dedrlx, dedrly, dedrlz, den, dhdl, dhdx, &
1330  dp1, dtau, dv2dz, dv2ijx, dv2ijy, dv2ijz, dv2j, dv3dz, dv3l, dv3ljx, dv3ljy, dv3ljz, &
1331  dv3lkx, dv3lky, dv3lkz, dv3rij, dv3rijx, dv3rijy, dv3rijz, dv3rik, dv3rikx, dv3riky, &
1332  dv3rikz, dwinv, dx, dxdz, dy, dz, ener_iat, fjx, fjy, fjz, fkx, fky, fkz, fz, h, lcos, &
1333  muhalf, par_a, par_alp, par_b, par_bet, par_bg, par_c, par_cap_a, par_cap_b, par_delta, &
1334  par_eta, par_gam, par_lam, par_mu, par_palp, par_qo, par_rh, par_sig, pz, qort, r, rinv, &
1335  rmainv, rmbinv, tau, temp0, temp1, u1, u2, u3, u4, u5, winv, x, xarg
1336  REAL(kind=dp) :: xinv, xinv3, z
1337 
1338 ! size of s2[]
1339 ! atom ID numbers for s2[]
1340 ! size of s3[]
1341 ! atom ID numbers for s3[]
1342 ! size of sz[]
1343 ! atom ID numbers for sz[]
1344 ! indices for the store arrays
1345 ! EDIP parameters
1346 
1347  par_cap_a = 5.6714030e0_dp
1348  par_cap_b = 2.0002804e0_dp
1349  par_rh = 1.2085196e0_dp
1350  par_a = 3.1213820e0_dp
1351  par_sig = 0.5774108e0_dp
1352  par_lam = 1.4533108e0_dp
1353  par_gam = 1.1247945e0_dp
1354  par_b = 3.1213820e0_dp
1355  par_c = 2.5609104e0_dp
1356  par_delta = 78.7590539e0_dp
1357  par_mu = 0.6966326e0_dp
1358  par_qo = 312.1341346e0_dp
1359  par_palp = 1.4074424e0_dp
1360  par_bet = 0.0070975e0_dp
1361  par_alp = 3.1083847e0_dp
1362 
1363  u1 = -0.165799e0_dp
1364  u2 = 32.557e0_dp
1365  u3 = 0.286198e0_dp
1366  u4 = 0.66e0_dp
1367 
1368  par_bg = par_a
1369  par_eta = par_delta/par_qo
1370 
1371  DO i = 1, nat
1372  ff(1, i) = 0.0e0_dp
1373  ff(2, i) = 0.0e0_dp
1374  ff(3, i) = 0.0e0_dp
1375  END DO
1376 
1377  coord = 0.e0_dp
1378  coord2 = 0.e0_dp
1379  ener = 0.e0_dp
1380  ener2 = 0.e0_dp
1381  istop = 0
1382 
1383 ! COMBINE COEFFICIENTS
1384 
1385  qort = sqrt(par_qo)
1386  muhalf = par_mu*0.5e0_dp
1387  u5 = u2*u4
1388  bmc = par_b - par_c
1389  cmbinv = 1.0e0_dp/(par_c - par_b)
1390 
1391 ! --- LEVEL 1: OUTER LOOP OVER ATOMS ---
1392 
1393  atoms: DO i = iat1, iat2
1394 
1395 ! RESET COORDINATION AND NEIGHBOR NUMBERS
1396 
1397  coord_iat = 0.e0_dp
1398  ener_iat = 0.e0_dp
1399  z = 0.0e0_dp
1400  n2 = 1
1401  n3 = 1
1402  nz = 1
1403 
1404 ! --- LEVEL 2: LOOP PREPASS OVER PAIRS ---
1405 
1406  DO n = lsta(1, i), lsta(2, i)
1407  j = lstb(n)
1408 
1409 ! PARTS OF TWO-BODY INTERACTION r<par_a
1410 
1411  num2(n2) = j
1412  dx = -rel(1, n)
1413  dy = -rel(2, n)
1414  dz = -rel(3, n)
1415  r = rel(4, n)
1416  rinv = rel(5, n)
1417  rmainv = 1.e0_dp/(r - par_a)
1418  s2_t0(n2) = par_cap_a*exp(par_sig*rmainv)
1419  s2_t1(n2) = (par_cap_b*rinv)**par_rh
1420  s2_t2(n2) = par_rh*rinv
1421  s2_t3(n2) = par_sig*rmainv*rmainv
1422  s2_dx(n2) = dx
1423  s2_dy(n2) = dy
1424  s2_dz(n2) = dz
1425  s2_r(n2) = r
1426  n2 = n2 + 1
1427  IF (n2 .GT. max_nbrs) THEN
1428  WRITE (*, *) 'WARNING enlarge max_nbrs'
1429  istop = 1
1430  RETURN
1431  END IF
1432 
1433 ! coordination number calculated with soft cutoff between first
1434 ! nearest neighbor and midpoint of first and second nearest neighbor
1435  IF (r .LE. 2.36e0_dp) THEN
1436  coord_iat = coord_iat + 1.e0_dp
1437  ELSE IF (r .GE. 3.12e0_dp) THEN
1438  ELSE
1439  xarg = (r - 2.36e0_dp)*(1.e0_dp/(3.12e0_dp - 2.36e0_dp))
1440  coord_iat = coord_iat + (2*xarg + 1.e0_dp)*(xarg - 1.e0_dp)**2
1441  END IF
1442 
1443 ! RADIAL PARTS OF THREE-BODY INTERACTION r<par_b
1444 
1445  IF (r .LT. par_bg) THEN
1446 
1447  num3(n3) = j
1448  rmbinv = 1.e0_dp/(r - par_bg)
1449  temp1 = par_gam*rmbinv
1450  temp0 = exp(temp1)
1451  s3_g(n3) = temp0
1452  s3_dg(n3) = -rmbinv*temp1*temp0
1453  s3_dx(n3) = dx
1454  s3_dy(n3) = dy
1455  s3_dz(n3) = dz
1456  s3_rinv(n3) = rinv
1457  s3_r(n3) = r
1458  n3 = n3 + 1
1459  IF (n3 .GT. max_nbrs) THEN
1460  WRITE (*, *) 'WARNING enlarge max_nbrs'
1461  istop = 1
1462  RETURN
1463  END IF
1464 
1465 ! COORDINATION AND NEIGHBOR FUNCTION par_c<r<par_b
1466 
1467  IF (r .LT. par_b) THEN
1468  IF (r .LT. par_c) THEN
1469  z = z + 1.e0_dp
1470  ELSE
1471  xinv = bmc/(r - par_c)
1472  xinv3 = xinv*xinv*xinv
1473  den = 1.e0_dp/(1 - xinv3)
1474  temp1 = par_alp*den
1475  fz = exp(temp1)
1476  z = z + fz
1477  numz(nz) = j
1478  sz_df(nz) = fz*temp1*den*3.e0_dp*xinv3*xinv*cmbinv
1479 ! df/dr
1480  sz_dx(nz) = dx
1481  sz_dy(nz) = dy
1482  sz_dz(nz) = dz
1483  sz_r(nz) = r
1484  nz = nz + 1
1485  IF (nz .GT. max_nbrs) THEN
1486  WRITE (*, *) 'WARNING enlarge max_nbrs'
1487  istop = 1
1488  RETURN
1489  END IF
1490  END IF
1491 ! r < par_C
1492  END IF
1493 ! r < par_b
1494  END IF
1495 ! r < par_bg
1496  END DO
1497 
1498 ! ZERO ACCUMULATION ARRAY FOR ENVIRONMENT FORCES
1499 
1500  DO nl = 1, nz - 1
1501  sz_sum(nl) = 0.e0_dp
1502  END DO
1503 
1504 ! ENVIRONMENT-DEPENDENCE OF PAIR INTERACTION
1505 
1506  temp0 = par_bet*z
1507  pz = par_palp*exp(-temp0*z)
1508 ! bond order
1509  dp1 = -2.e0_dp*temp0*pz
1510 ! derivative of bond order
1511 
1512 ! --- LEVEL 2: LOOP FOR PAIR INTERACTIONS ---
1513 
1514  DO nj = 1, n2 - 1
1515 
1516  temp0 = s2_t1(nj) - pz
1517 
1518 ! two-body energy V2(rij,Z)
1519 
1520  ener_iat = ener_iat + temp0*s2_t0(nj)
1521 
1522 ! two-body forces
1523 
1524  dv2j = -s2_t0(nj)*(s2_t1(nj)*s2_t2(nj) + temp0*s2_t3(nj))
1525 ! dV2/dr
1526  dv2ijx = dv2j*s2_dx(nj)
1527  dv2ijy = dv2j*s2_dy(nj)
1528  dv2ijz = dv2j*s2_dz(nj)
1529  ff(1, i) = ff(1, i) + dv2ijx
1530  ff(2, i) = ff(2, i) + dv2ijy
1531  ff(3, i) = ff(3, i) + dv2ijz
1532  j = num2(nj)
1533  ff(1, j) = ff(1, j) - dv2ijx
1534  ff(2, j) = ff(2, j) - dv2ijy
1535  ff(3, j) = ff(3, j) - dv2ijz
1536 
1537 ! --- LEVEL 3: LOOP FOR PAIR COORDINATION FORCES ---
1538 
1539  dv2dz = -dp1*s2_t0(nj)
1540  DO nl = 1, nz - 1
1541  sz_sum(nl) = sz_sum(nl) + dv2dz
1542  END DO
1543 
1544  END DO
1545 
1546 ! COORDINATION-DEPENDENCE OF THREE-BODY INTERACTION
1547 
1548  winv = qort*exp(-muhalf*z)
1549 ! inverse width of angular function
1550  dwinv = -muhalf*winv
1551 ! its derivative
1552  temp0 = exp(-u4*z)
1553  tau = u1 + u2*temp0*(u3 - temp0)
1554 ! -cosine of angular minimum
1555  dtau = u5*temp0*(2*temp0 - u3)
1556 ! its derivative
1557 
1558 ! --- LEVEL 2: FIRST LOOP FOR THREE-BODY INTERACTIONS ---
1559 
1560  DO nj = 1, n3 - 2
1561 
1562  j = num3(nj)
1563 
1564 ! --- LEVEL 3: SECOND LOOP FOR THREE-BODY INTERACTIONS ---
1565 
1566  DO nk = nj + 1, n3 - 1
1567 
1568  k = num3(nk)
1569 
1570 ! angular function h(l,Z)
1571 
1572  lcos = s3_dx(nj)*s3_dx(nk) + s3_dy(nj)*s3_dy(nk) + s3_dz(nj)*s3_dz(nk)
1573  x = (lcos + tau)*winv
1574  temp0 = exp(-x*x)
1575 
1576  h = par_lam*(1 - temp0 + par_eta*x*x)
1577  dhdx = 2*par_lam*x*(temp0 + par_eta)
1578 
1579  dhdl = dhdx*winv
1580 
1581 ! three-body energy
1582 
1583  temp1 = s3_g(nj)*s3_g(nk)
1584  ener_iat = ener_iat + temp1*h
1585 
1586 ! (-) radial force on atom j
1587 
1588  dv3rij = s3_dg(nj)*s3_g(nk)*h
1589  dv3rijx = dv3rij*s3_dx(nj)
1590  dv3rijy = dv3rij*s3_dy(nj)
1591  dv3rijz = dv3rij*s3_dz(nj)
1592  fjx = dv3rijx
1593  fjy = dv3rijy
1594  fjz = dv3rijz
1595 
1596 ! (-) radial force on atom k
1597 
1598  dv3rik = s3_g(nj)*s3_dg(nk)*h
1599  dv3rikx = dv3rik*s3_dx(nk)
1600  dv3riky = dv3rik*s3_dy(nk)
1601  dv3rikz = dv3rik*s3_dz(nk)
1602  fkx = dv3rikx
1603  fky = dv3riky
1604  fkz = dv3rikz
1605 
1606 ! (-) angular force on j
1607 
1608  dv3l = temp1*dhdl
1609  dv3ljx = dv3l*(s3_dx(nk) - lcos*s3_dx(nj))*s3_rinv(nj)
1610  dv3ljy = dv3l*(s3_dy(nk) - lcos*s3_dy(nj))*s3_rinv(nj)
1611  dv3ljz = dv3l*(s3_dz(nk) - lcos*s3_dz(nj))*s3_rinv(nj)
1612  fjx = fjx + dv3ljx
1613  fjy = fjy + dv3ljy
1614  fjz = fjz + dv3ljz
1615 
1616 ! (-) angular force on k
1617 
1618  dv3lkx = dv3l*(s3_dx(nj) - lcos*s3_dx(nk))*s3_rinv(nk)
1619  dv3lky = dv3l*(s3_dy(nj) - lcos*s3_dy(nk))*s3_rinv(nk)
1620  dv3lkz = dv3l*(s3_dz(nj) - lcos*s3_dz(nk))*s3_rinv(nk)
1621  fkx = fkx + dv3lkx
1622  fky = fky + dv3lky
1623  fkz = fkz + dv3lkz
1624 
1625 ! apply radial + angular forces to i, j, k
1626 
1627  ff(1, j) = ff(1, j) - fjx
1628  ff(2, j) = ff(2, j) - fjy
1629  ff(3, j) = ff(3, j) - fjz
1630  ff(1, k) = ff(1, k) - fkx
1631  ff(2, k) = ff(2, k) - fky
1632  ff(3, k) = ff(3, k) - fkz
1633  ff(1, i) = ff(1, i) + fjx + fkx
1634  ff(2, i) = ff(2, i) + fjy + fky
1635  ff(3, i) = ff(3, i) + fjz + fkz
1636 
1637 ! prefactor for 4-body forces from coordination
1638  dxdz = dwinv*(lcos + tau) + winv*dtau
1639  dv3dz = temp1*dhdx*dxdz
1640 
1641 ! --- LEVEL 4: LOOP FOR THREE-BODY COORDINATION FORCES ---
1642 
1643  DO nl = 1, nz - 1
1644  sz_sum(nl) = sz_sum(nl) + dv3dz
1645  END DO
1646  END DO
1647  END DO
1648 
1649 ! --- LEVEL 2: LOOP TO APPLY COORDINATION FORCES ---
1650 
1651  DO nl = 1, nz - 1
1652 
1653  dedrl = sz_sum(nl)*sz_df(nl)
1654  dedrlx = dedrl*sz_dx(nl)
1655  dedrly = dedrl*sz_dy(nl)
1656  dedrlz = dedrl*sz_dz(nl)
1657  ff(1, i) = ff(1, i) + dedrlx
1658  ff(2, i) = ff(2, i) + dedrly
1659  ff(3, i) = ff(3, i) + dedrlz
1660  l = numz(nl)
1661  ff(1, l) = ff(1, l) - dedrlx
1662  ff(2, l) = ff(2, l) - dedrly
1663  ff(3, l) = ff(3, l) - dedrlz
1664 
1665  END DO
1666 
1667  coord = coord + coord_iat
1668  coord2 = coord2 + coord_iat**2
1669  ener = ener + ener_iat
1670  ener2 = ener2 + ener_iat**2
1671 
1672  END DO atoms
1673 
1674  RETURN
1675  END SUBROUTINE subfeniat_b
1676 
1677 ! **************************************************************************************************
1678 !> \brief ...
1679 !> \param iat ...
1680 !> \param nn ...
1681 !> \param ncx ...
1682 !> \param ll1 ...
1683 !> \param ll2 ...
1684 !> \param ll3 ...
1685 !> \param l1 ...
1686 !> \param l2 ...
1687 !> \param l3 ...
1688 !> \param myspace ...
1689 !> \param rxyz ...
1690 !> \param icell ...
1691 !> \param lstb ...
1692 !> \param lay ...
1693 !> \param rel ...
1694 !> \param cut2 ...
1695 !> \param indlst ...
1696 ! **************************************************************************************************
1697  SUBROUTINE sublstiat_b(iat, nn, ncx, ll1, ll2, ll3, l1, l2, l3, myspace, &
1698  rxyz, icell, lstb, lay, rel, cut2, indlst)
1699 ! finds the neighbours of atom iat (specified by lsta and lstb) and and
1700 ! the relative position rel of iat with respect to these neighbours
1701  INTEGER :: iat, nn, ncx, ll1, ll2, ll3, l1, l2, l3, &
1702  myspace
1703  REAL(kind=dp) :: rxyz
1704  INTEGER :: icell, lstb, lay
1705  REAL(kind=dp) :: rel, cut2
1706  INTEGER :: indlst
1707 
1708  dimension rxyz(3, nn), lay(nn), icell(0:ncx, -1:ll1, -1:ll2, -1:ll3), &
1709  lstb(0:myspace - 1), rel(5, 0:myspace - 1)
1710 
1711  INTEGER :: jat, k1, k2, k3, jj
1712  REAL(kind=dp) :: rr2, tt, tti, xrel, yrel, zrel
1713 
1714  DO k3 = l3 - 1, l3 + 1
1715  DO k2 = l2 - 1, l2 + 1
1716  DO k1 = l1 - 1, l1 + 1
1717  DO jj = 1, icell(0, k1, k2, k3)
1718  jat = icell(jj, k1, k2, k3)
1719  IF (jat .EQ. iat) cycle
1720  xrel = rxyz(1, iat) - rxyz(1, jat)
1721  yrel = rxyz(2, iat) - rxyz(2, jat)
1722  zrel = rxyz(3, iat) - rxyz(3, jat)
1723  rr2 = xrel**2 + yrel**2 + zrel**2
1724  IF (rr2 .LE. cut2) THEN
1725  indlst = min(indlst, myspace - 1)
1726  lstb(indlst) = lay(jat)
1727 ! write(*,*) 'iat,indlst,lay(jat)',iat,indlst,lay(jat)
1728  tt = sqrt(rr2)
1729  tti = 1.e0_dp/tt
1730  rel(1, indlst) = xrel*tti
1731  rel(2, indlst) = yrel*tti
1732  rel(3, indlst) = zrel*tti
1733  rel(4, indlst) = tt
1734  rel(5, indlst) = tti
1735  indlst = indlst + 1
1736  END IF
1737  END DO
1738  END DO
1739  END DO
1740  END DO
1741 
1742  RETURN
1743  END SUBROUTINE sublstiat_b
1744 
1745 ! **************************************************************************************************
1746 !> \brief Lenosky's "highly optimized empirical potential model of silicon"
1747 !> by Stefan Goedecker
1748 !> \param nat number of atoms
1749 !> \param alat lattice constants of the orthorombic box containing the particles
1750 !> \param rxyz0 atomic positions in Angstrom, may be modified on output.
1751 !> If an atom is outside the box the program will bring it back
1752 !> into the box by translations through alat
1753 !> \param fxyz forces in eV/A
1754 !> \param ener total energy in eV
1755 !> \param coord average coordination number
1756 !> \param ener_var variance of the energy/atom
1757 !> \param coord_var variance of the coordination number
1758 !> \param count count is increased by one per call, has to be initialized
1759 !> to 0.e0_dp before first call of eip_bazant
1760 !> \par Literature
1761 !> T. Lenosky, et. al.: Highly optimized empirical potential model of silicon;
1762 !> Modeling Simul. Sci. Eng., 8 (2000)
1763 !> S. Goedecker: Optimization and parallelization of a force field for silicon
1764 !> using OpenMP; CPC 148, 1 (2002)
1765 !> \par History
1766 !> 03.2006 initial create [tdk]
1767 !> \author Thomas D. Kuehne (tkuehne@phys.chem.ethz.ch)
1768 ! **************************************************************************************************
1769  SUBROUTINE eip_lenosky_silicon(nat, alat, rxyz0, fxyz, ener, coord, ener_var, &
1770  coord_var, count)
1771 
1772  INTEGER :: nat
1773  REAL(kind=dp) :: alat, rxyz0, fxyz, ener, coord, &
1774  ener_var, coord_var, count
1775 
1776  dimension rxyz0(3, nat), fxyz(3, nat), alat(3)
1777  REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :) :: rxyz
1778  INTEGER, ALLOCATABLE, DIMENSION(:, :) :: lsta
1779  INTEGER, ALLOCATABLE, DIMENSION(:) :: lstb
1780  INTEGER, ALLOCATABLE, DIMENSION(:) :: lay
1781  INTEGER, ALLOCATABLE, DIMENSION(:, :, :, :) :: icell
1782  REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :) :: rel
1783  REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :) :: txyz
1784  REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :) :: f2ij, f3ij, f3ik
1785 
1786  REAL(kind=dp) :: coord2, cut, cut2, ener2, tcoord, &
1787  tcoord2, tener, tener2
1788  INTEGER :: i, iam, iat, iat1, iat2, ii, il, in, indlst, indlstx, &
1789  istop, istopg, l1, l2, l3, ll1, ll2, ll3, lot, ncx, nn, &
1790  nnbrx, npjkx, npjx, laymx, npr, rlc1i, rlc2i, rlc3i, &
1791  myspace, myspaceout
1792 
1793 ! tmax_phi= 0.4500000e+01_dp
1794 ! cut=tmax_phi
1795  cut = 0.4500000e+01_dp
1796 
1797  IF (count .EQ. 0) OPEN (unit=10, file='lenosky.mon', status='unknown')
1798  count = count + 1.e0_dp
1799 
1800 ! linear scaling calculation of verlet list
1801  ll1 = int(alat(1)/cut)
1802  IF (ll1 .LT. 1) cpabort("alat(1) too small")
1803  ll2 = int(alat(2)/cut)
1804  IF (ll2 .LT. 1) cpabort("alat(2) too small")
1805  ll3 = int(alat(3)/cut)
1806  IF (ll3 .LT. 1) cpabort("alat(3) too small")
1807 
1808 ! determine number of threads
1809  npr = 1
1810 !$OMP PARALLEL PRIVATE(iam) SHARED (npr) DEFAULT(NONE)
1811 !$ iam = omp_get_thread_num()
1812 !$ if (iam .eq. 0) npr = omp_get_num_threads()
1813 !$OMP END PARALLEL
1814 
1815 ! linear scaling calculation of verlet list
1816 
1817  IF (npr .LE. 1) THEN !serial if too few processors to gain by parallelizing
1818 
1819 ! set ncx for serial case, ncx for parallel case set below
1820  ncx = 16
1821  loop_ncx_s: DO
1822  ALLOCATE (icell(0:ncx, -1:ll1, -1:ll2, -1:ll3))
1823  icell(0, -1:ll1, -1:ll2, -1:ll3) = 0
1824  rlc1i = int(ll1/alat(1))
1825  rlc2i = int(ll2/alat(2))
1826  rlc3i = int(ll3/alat(3))
1827 
1828  loop_iat_s: DO iat = 1, nat
1829  rxyz0(1, iat) = modulo(modulo(rxyz0(1, iat), alat(1)), alat(1))
1830  rxyz0(2, iat) = modulo(modulo(rxyz0(2, iat), alat(2)), alat(2))
1831  rxyz0(3, iat) = modulo(modulo(rxyz0(3, iat), alat(3)), alat(3))
1832  l1 = int(rxyz0(1, iat)*rlc1i)
1833  l2 = int(rxyz0(2, iat)*rlc2i)
1834  l3 = int(rxyz0(3, iat)*rlc3i)
1835 
1836  ii = icell(0, l1, l2, l3)
1837  ii = ii + 1
1838  icell(0, l1, l2, l3) = ii
1839  IF (ii .GT. ncx) THEN
1840  WRITE (10, *) count, 'NCX too small', ncx
1841  DEALLOCATE (icell)
1842  ncx = ncx*2
1843  cycle loop_ncx_s
1844  END IF
1845  icell(ii, l1, l2, l3) = iat
1846  END DO loop_iat_s
1847  EXIT loop_ncx_s
1848  END DO loop_ncx_s
1849 
1850  ELSE ! parallel case
1851 
1852 ! periodization of particles can be done in parallel
1853 !$OMP PARALLEL DO SHARED (alat,nat,rxyz0) PRIVATE(iat) DEFAULT(NONE)
1854  DO iat = 1, nat
1855  rxyz0(1, iat) = modulo(modulo(rxyz0(1, iat), alat(1)), alat(1))
1856  rxyz0(2, iat) = modulo(modulo(rxyz0(2, iat), alat(2)), alat(2))
1857  rxyz0(3, iat) = modulo(modulo(rxyz0(3, iat), alat(3)), alat(3))
1858  END DO
1859 !$OMP END PARALLEL DO
1860 
1861 ! assignment to cell is done serially
1862 ! set ncx for parallel case, ncx for serial case set above
1863  ncx = 16
1864  loop_ncx_p: DO
1865  ALLOCATE (icell(0:ncx, -1:ll1, -1:ll2, -1:ll3))
1866  icell(0, -1:ll1, -1:ll2, -1:ll3) = 0
1867  rlc1i = int(ll1/alat(1))
1868  rlc2i = int(ll2/alat(2))
1869  rlc3i = int(ll3/alat(3))
1870 
1871  loop_iat_p: DO iat = 1, nat
1872  l1 = int(rxyz0(1, iat)*rlc1i)
1873  l2 = int(rxyz0(2, iat)*rlc2i)
1874  l3 = int(rxyz0(3, iat)*rlc3i)
1875  ii = icell(0, l1, l2, l3)
1876  ii = ii + 1
1877  icell(0, l1, l2, l3) = ii
1878  IF (ii .GT. ncx) THEN
1879  WRITE (10, *) count, 'NCX too small', ncx
1880  DEALLOCATE (icell)
1881  ncx = ncx*2
1882  cycle loop_ncx_p
1883  END IF
1884  icell(ii, l1, l2, l3) = iat
1885  END DO loop_iat_p
1886  EXIT loop_ncx_p
1887  END DO loop_ncx_p
1888 
1889  END IF
1890 
1891 ! duplicate all atoms within boundary layer
1892  laymx = ncx*(2*ll1*ll2 + 2*ll1*ll3 + 2*ll2*ll3 + 4*ll1 + 4*ll2 + 4*ll3 + 8)
1893  nn = nat + laymx
1894  ALLOCATE (rxyz(3, nn), lay(nn))
1895  DO iat = 1, nat
1896  lay(iat) = iat
1897  rxyz(1, iat) = rxyz0(1, iat)
1898  rxyz(2, iat) = rxyz0(2, iat)
1899  rxyz(3, iat) = rxyz0(3, iat)
1900  END DO
1901  il = nat
1902 ! xy plane
1903  DO l2 = 0, ll2 - 1
1904  DO l1 = 0, ll1 - 1
1905 
1906  in = icell(0, l1, l2, 0)
1907  icell(0, l1, l2, ll3) = in
1908  DO ii = 1, in
1909  i = icell(ii, l1, l2, 0)
1910  il = il + 1
1911  IF (il .GT. nn) cpabort("enlarge laymx")
1912  lay(il) = i
1913  icell(ii, l1, l2, ll3) = il
1914  rxyz(1, il) = rxyz(1, i)
1915  rxyz(2, il) = rxyz(2, i)
1916  rxyz(3, il) = rxyz(3, i) + alat(3)
1917  END DO
1918 
1919  in = icell(0, l1, l2, ll3 - 1)
1920  icell(0, l1, l2, -1) = in
1921  DO ii = 1, in
1922  i = icell(ii, l1, l2, ll3 - 1)
1923  il = il + 1
1924  IF (il .GT. nn) cpabort("enlarge laymx")
1925  lay(il) = i
1926  icell(ii, l1, l2, -1) = il
1927  rxyz(1, il) = rxyz(1, i)
1928  rxyz(2, il) = rxyz(2, i)
1929  rxyz(3, il) = rxyz(3, i) - alat(3)
1930  END DO
1931 
1932  END DO
1933  END DO
1934 
1935 ! yz plane
1936  DO l3 = 0, ll3 - 1
1937  DO l2 = 0, ll2 - 1
1938 
1939  in = icell(0, 0, l2, l3)
1940  icell(0, ll1, l2, l3) = in
1941  DO ii = 1, in
1942  i = icell(ii, 0, l2, l3)
1943  il = il + 1
1944  IF (il .GT. nn) cpabort("enlarge laymx")
1945  lay(il) = i
1946  icell(ii, ll1, l2, l3) = il
1947  rxyz(1, il) = rxyz(1, i) + alat(1)
1948  rxyz(2, il) = rxyz(2, i)
1949  rxyz(3, il) = rxyz(3, i)
1950  END DO
1951 
1952  in = icell(0, ll1 - 1, l2, l3)
1953  icell(0, -1, l2, l3) = in
1954  DO ii = 1, in
1955  i = icell(ii, ll1 - 1, l2, l3)
1956  il = il + 1
1957  IF (il .GT. nn) cpabort("enlarge laymx")
1958  lay(il) = i
1959  icell(ii, -1, l2, l3) = il
1960  rxyz(1, il) = rxyz(1, i) - alat(1)
1961  rxyz(2, il) = rxyz(2, i)
1962  rxyz(3, il) = rxyz(3, i)
1963  END DO
1964 
1965  END DO
1966  END DO
1967 
1968 ! xz plane
1969  DO l3 = 0, ll3 - 1
1970  DO l1 = 0, ll1 - 1
1971 
1972  in = icell(0, l1, 0, l3)
1973  icell(0, l1, ll2, l3) = in
1974  DO ii = 1, in
1975  i = icell(ii, l1, 0, l3)
1976  il = il + 1
1977  IF (il .GT. nn) cpabort("enlarge laymx")
1978  lay(il) = i
1979  icell(ii, l1, ll2, l3) = il
1980  rxyz(1, il) = rxyz(1, i)
1981  rxyz(2, il) = rxyz(2, i) + alat(2)
1982  rxyz(3, il) = rxyz(3, i)
1983  END DO
1984 
1985  in = icell(0, l1, ll2 - 1, l3)
1986  icell(0, l1, -1, l3) = in
1987  DO ii = 1, in
1988  i = icell(ii, l1, ll2 - 1, l3)
1989  il = il + 1
1990  IF (il .GT. nn) cpabort("enlarge laymx")
1991  lay(il) = i
1992  icell(ii, l1, -1, l3) = il
1993  rxyz(1, il) = rxyz(1, i)
1994  rxyz(2, il) = rxyz(2, i) - alat(2)
1995  rxyz(3, il) = rxyz(3, i)
1996  END DO
1997 
1998  END DO
1999  END DO
2000 
2001 ! x axis
2002  DO l1 = 0, ll1 - 1
2003 
2004  in = icell(0, l1, 0, 0)
2005  icell(0, l1, ll2, ll3) = in
2006  DO ii = 1, in
2007  i = icell(ii, l1, 0, 0)
2008  il = il + 1
2009  IF (il .GT. nn) cpabort("enlarge laymx")
2010  lay(il) = i
2011  icell(ii, l1, ll2, ll3) = il
2012  rxyz(1, il) = rxyz(1, i)
2013  rxyz(2, il) = rxyz(2, i) + alat(2)
2014  rxyz(3, il) = rxyz(3, i) + alat(3)
2015  END DO
2016 
2017  in = icell(0, l1, 0, ll3 - 1)
2018  icell(0, l1, ll2, -1) = in
2019  DO ii = 1, in
2020  i = icell(ii, l1, 0, ll3 - 1)
2021  il = il + 1
2022  IF (il .GT. nn) cpabort("enlarge laymx")
2023  lay(il) = i
2024  icell(ii, l1, ll2, -1) = il
2025  rxyz(1, il) = rxyz(1, i)
2026  rxyz(2, il) = rxyz(2, i) + alat(2)
2027  rxyz(3, il) = rxyz(3, i) - alat(3)
2028  END DO
2029 
2030  in = icell(0, l1, ll2 - 1, 0)
2031  icell(0, l1, -1, ll3) = in
2032  DO ii = 1, in
2033  i = icell(ii, l1, ll2 - 1, 0)
2034  il = il + 1
2035  IF (il .GT. nn) cpabort("enlarge laymx")
2036  lay(il) = i
2037  icell(ii, l1, -1, ll3) = il
2038  rxyz(1, il) = rxyz(1, i)
2039  rxyz(2, il) = rxyz(2, i) - alat(2)
2040  rxyz(3, il) = rxyz(3, i) + alat(3)
2041  END DO
2042 
2043  in = icell(0, l1, ll2 - 1, ll3 - 1)
2044  icell(0, l1, -1, -1) = in
2045  DO ii = 1, in
2046  i = icell(ii, l1, ll2 - 1, ll3 - 1)
2047  il = il + 1
2048  IF (il .GT. nn) cpabort("enlarge laymx")
2049  lay(il) = i
2050  icell(ii, l1, -1, -1) = il
2051  rxyz(1, il) = rxyz(1, i)
2052  rxyz(2, il) = rxyz(2, i) - alat(2)
2053  rxyz(3, il) = rxyz(3, i) - alat(3)
2054  END DO
2055 
2056  END DO
2057 
2058 ! y axis
2059  DO l2 = 0, ll2 - 1
2060 
2061  in = icell(0, 0, l2, 0)
2062  icell(0, ll1, l2, ll3) = in
2063  DO ii = 1, in
2064  i = icell(ii, 0, l2, 0)
2065  il = il + 1
2066  IF (il .GT. nn) cpabort("enlarge laymx")
2067  lay(il) = i
2068  icell(ii, ll1, l2, ll3) = il
2069  rxyz(1, il) = rxyz(1, i) + alat(1)
2070  rxyz(2, il) = rxyz(2, i)
2071  rxyz(3, il) = rxyz(3, i) + alat(3)
2072  END DO
2073 
2074  in = icell(0, 0, l2, ll3 - 1)
2075  icell(0, ll1, l2, -1) = in
2076  DO ii = 1, in
2077  i = icell(ii, 0, l2, ll3 - 1)
2078  il = il + 1
2079  IF (il .GT. nn) cpabort("enlarge laymx")
2080  lay(il) = i
2081  icell(ii, ll1, l2, -1) = il
2082  rxyz(1, il) = rxyz(1, i) + alat(1)
2083  rxyz(2, il) = rxyz(2, i)
2084  rxyz(3, il) = rxyz(3, i) - alat(3)
2085  END DO
2086 
2087  in = icell(0, ll1 - 1, l2, 0)
2088  icell(0, -1, l2, ll3) = in
2089  DO ii = 1, in
2090  i = icell(ii, ll1 - 1, l2, 0)
2091  il = il + 1
2092  IF (il .GT. nn) cpabort("enlarge laymx")
2093  lay(il) = i
2094  icell(ii, -1, l2, ll3) = il
2095  rxyz(1, il) = rxyz(1, i) - alat(1)
2096  rxyz(2, il) = rxyz(2, i)
2097  rxyz(3, il) = rxyz(3, i) + alat(3)
2098  END DO
2099 
2100  in = icell(0, ll1 - 1, l2, ll3 - 1)
2101  icell(0, -1, l2, -1) = in
2102  DO ii = 1, in
2103  i = icell(ii, ll1 - 1, l2, ll3 - 1)
2104  il = il + 1
2105  IF (il .GT. nn) cpabort("enlarge laymx")
2106  lay(il) = i
2107  icell(ii, -1, l2, -1) = il
2108  rxyz(1, il) = rxyz(1, i) - alat(1)
2109  rxyz(2, il) = rxyz(2, i)
2110  rxyz(3, il) = rxyz(3, i) - alat(3)
2111  END DO
2112 
2113  END DO
2114 
2115 ! z axis
2116  DO l3 = 0, ll3 - 1
2117 
2118  in = icell(0, 0, 0, l3)
2119  icell(0, ll1, ll2, l3) = in
2120  DO ii = 1, in
2121  i = icell(ii, 0, 0, l3)
2122  il = il + 1
2123  IF (il .GT. nn) cpabort("enlarge laymx")
2124  lay(il) = i
2125  icell(ii, ll1, ll2, l3) = il
2126  rxyz(1, il) = rxyz(1, i) + alat(1)
2127  rxyz(2, il) = rxyz(2, i) + alat(2)
2128  rxyz(3, il) = rxyz(3, i)
2129  END DO
2130 
2131  in = icell(0, ll1 - 1, 0, l3)
2132  icell(0, -1, ll2, l3) = in
2133  DO ii = 1, in
2134  i = icell(ii, ll1 - 1, 0, l3)
2135  il = il + 1
2136  IF (il .GT. nn) cpabort("enlarge laymx")
2137  lay(il) = i
2138  icell(ii, -1, ll2, l3) = il
2139  rxyz(1, il) = rxyz(1, i) - alat(1)
2140  rxyz(2, il) = rxyz(2, i) + alat(2)
2141  rxyz(3, il) = rxyz(3, i)
2142  END DO
2143 
2144  in = icell(0, 0, ll2 - 1, l3)
2145  icell(0, ll1, -1, l3) = in
2146  DO ii = 1, in
2147  i = icell(ii, 0, ll2 - 1, l3)
2148  il = il + 1
2149  IF (il .GT. nn) cpabort("enlarge laymx")
2150  lay(il) = i
2151  icell(ii, ll1, -1, l3) = il
2152  rxyz(1, il) = rxyz(1, i) + alat(1)
2153  rxyz(2, il) = rxyz(2, i) - alat(2)
2154  rxyz(3, il) = rxyz(3, i)
2155  END DO
2156 
2157  in = icell(0, ll1 - 1, ll2 - 1, l3)
2158  icell(0, -1, -1, l3) = in
2159  DO ii = 1, in
2160  i = icell(ii, ll1 - 1, ll2 - 1, l3)
2161  il = il + 1
2162  IF (il .GT. nn) cpabort("enlarge laymx")
2163  lay(il) = i
2164  icell(ii, -1, -1, l3) = il
2165  rxyz(1, il) = rxyz(1, i) - alat(1)
2166  rxyz(2, il) = rxyz(2, i) - alat(2)
2167  rxyz(3, il) = rxyz(3, i)
2168  END DO
2169 
2170  END DO
2171 
2172 ! corners
2173  in = icell(0, 0, 0, 0)
2174  icell(0, ll1, ll2, ll3) = in
2175  DO ii = 1, in
2176  i = icell(ii, 0, 0, 0)
2177  il = il + 1
2178  IF (il .GT. nn) cpabort("enlarge laymx")
2179  lay(il) = i
2180  icell(ii, ll1, ll2, ll3) = il
2181  rxyz(1, il) = rxyz(1, i) + alat(1)
2182  rxyz(2, il) = rxyz(2, i) + alat(2)
2183  rxyz(3, il) = rxyz(3, i) + alat(3)
2184  END DO
2185 
2186  in = icell(0, ll1 - 1, 0, 0)
2187  icell(0, -1, ll2, ll3) = in
2188  DO ii = 1, in
2189  i = icell(ii, ll1 - 1, 0, 0)
2190  il = il + 1
2191  IF (il .GT. nn) cpabort("enlarge laymx")
2192  lay(il) = i
2193  icell(ii, -1, ll2, ll3) = il
2194  rxyz(1, il) = rxyz(1, i) - alat(1)
2195  rxyz(2, il) = rxyz(2, i) + alat(2)
2196  rxyz(3, il) = rxyz(3, i) + alat(3)
2197  END DO
2198 
2199  in = icell(0, 0, ll2 - 1, 0)
2200  icell(0, ll1, -1, ll3) = in
2201  DO ii = 1, in
2202  i = icell(ii, 0, ll2 - 1, 0)
2203  il = il + 1
2204  IF (il .GT. nn) cpabort("enlarge laymx")
2205  lay(il) = i
2206  icell(ii, ll1, -1, ll3) = il
2207  rxyz(1, il) = rxyz(1, i) + alat(1)
2208  rxyz(2, il) = rxyz(2, i) - alat(2)
2209  rxyz(3, il) = rxyz(3, i) + alat(3)
2210  END DO
2211 
2212  in = icell(0, ll1 - 1, ll2 - 1, 0)
2213  icell(0, -1, -1, ll3) = in
2214  DO ii = 1, in
2215  i = icell(ii, ll1 - 1, ll2 - 1, 0)
2216  il = il + 1
2217  IF (il .GT. nn) cpabort("enlarge laymx")
2218  lay(il) = i
2219  icell(ii, -1, -1, ll3) = il
2220  rxyz(1, il) = rxyz(1, i) - alat(1)
2221  rxyz(2, il) = rxyz(2, i) - alat(2)
2222  rxyz(3, il) = rxyz(3, i) + alat(3)
2223  END DO
2224 
2225  in = icell(0, 0, 0, ll3 - 1)
2226  icell(0, ll1, ll2, -1) = in
2227  DO ii = 1, in
2228  i = icell(ii, 0, 0, ll3 - 1)
2229  il = il + 1
2230  IF (il .GT. nn) cpabort("enlarge laymx")
2231  lay(il) = i
2232  icell(ii, ll1, ll2, -1) = il
2233  rxyz(1, il) = rxyz(1, i) + alat(1)
2234  rxyz(2, il) = rxyz(2, i) + alat(2)
2235  rxyz(3, il) = rxyz(3, i) - alat(3)
2236  END DO
2237 
2238  in = icell(0, ll1 - 1, 0, ll3 - 1)
2239  icell(0, -1, ll2, -1) = in
2240  DO ii = 1, in
2241  i = icell(ii, ll1 - 1, 0, ll3 - 1)
2242  il = il + 1
2243  IF (il .GT. nn) cpabort("enlarge laymx")
2244  lay(il) = i
2245  icell(ii, -1, ll2, -1) = il
2246  rxyz(1, il) = rxyz(1, i) - alat(1)
2247  rxyz(2, il) = rxyz(2, i) + alat(2)
2248  rxyz(3, il) = rxyz(3, i) - alat(3)
2249  END DO
2250 
2251  in = icell(0, 0, ll2 - 1, ll3 - 1)
2252  icell(0, ll1, -1, -1) = in
2253  DO ii = 1, in
2254  i = icell(ii, 0, ll2 - 1, ll3 - 1)
2255  il = il + 1
2256  IF (il .GT. nn) cpabort("enlarge laymx")
2257  lay(il) = i
2258  icell(ii, ll1, -1, -1) = il
2259  rxyz(1, il) = rxyz(1, i) + alat(1)
2260  rxyz(2, il) = rxyz(2, i) - alat(2)
2261  rxyz(3, il) = rxyz(3, i) - alat(3)
2262  END DO
2263 
2264  in = icell(0, ll1 - 1, ll2 - 1, ll3 - 1)
2265  icell(0, -1, -1, -1) = in
2266  DO ii = 1, in
2267  i = icell(ii, ll1 - 1, ll2 - 1, ll3 - 1)
2268  il = il + 1
2269  IF (il .GT. nn) cpabort("enlarge laymx")
2270  lay(il) = i
2271  icell(ii, -1, -1, -1) = il
2272  rxyz(1, il) = rxyz(1, i) - alat(1)
2273  rxyz(2, il) = rxyz(2, i) - alat(2)
2274  rxyz(3, il) = rxyz(3, i) - alat(3)
2275  END DO
2276 
2277  ALLOCATE (lsta(2, nat))
2278  nnbrx = 36
2279  loop_nnbrx: DO
2280  ALLOCATE (lstb(nnbrx*nat), rel(5, nnbrx*nat))
2281 
2282  indlstx = 0
2283 
2284 !$OMP PARALLEL DEFAULT(NONE) &
2285 !$OMP PRIVATE(iat,cut2,iam,ii,indlst,l1,l2,l3,myspace,npr) &
2286 !$OMP SHARED (indlstx,nat,nn,nnbrx,ncx,ll1,ll2,ll3,icell,lsta,lstb,lay, &
2287 !$OMP rel,rxyz,cut,myspaceout)
2288 
2289  npr = 1
2290 !$ npr = omp_get_num_threads()
2291  iam = 0
2292 !$ iam = omp_get_thread_num()
2293 
2294  cut2 = cut**2
2295 ! assign contiguous portions of the arrays lstb and rel to the threads
2296  myspace = (nat*nnbrx)/npr
2297  IF (iam .EQ. 0) myspaceout = myspace
2298 ! Verlet list, relative positions
2299  indlst = 0
2300  loop_l3: DO l3 = 0, ll3 - 1
2301  loop_l2: DO l2 = 0, ll2 - 1
2302  loop_l1: DO l1 = 0, ll1 - 1
2303  loop_ii: DO ii = 1, icell(0, l1, l2, l3)
2304  iat = icell(ii, l1, l2, l3)
2305  IF (((iat - 1)*npr)/nat .EQ. iam) THEN
2306 ! write(*,*) 'sublstiat:iam,iat',iam,iat
2307  lsta(1, iat) = iam*myspace + indlst + 1
2308  CALL sublstiat_l(iat, nn, ncx, ll1, ll2, ll3, l1, l2, l3, myspace, &
2309  rxyz, icell, lstb(iam*myspace + 1), lay, &
2310  rel(1, iam*myspace + 1), cut2, indlst)
2311  lsta(2, iat) = iam*myspace + indlst
2312 ! write(*,'(a,4(x,i3),100(x,i2))') &
2313 ! 'iam,iat,lsta',iam,iat,lsta(1,iat),lsta(2,iat), &
2314 ! (lstb(j),j=lsta(1,iat),lsta(2,iat))
2315  END IF
2316  END DO loop_ii
2317  END DO loop_l1
2318  END DO loop_l2
2319  END DO loop_l3
2320 
2321 !$OMP CRITICAL
2322  indlstx = max(indlstx, indlst)
2323 !$OMP END CRITICAL
2324 !$OMP END PARALLEL
2325 
2326  IF (indlstx .GE. myspaceout) THEN
2327  WRITE (10, *) count, 'NNBRX too small', nnbrx
2328  DEALLOCATE (lstb, rel)
2329  nnbrx = 3*nnbrx/2
2330  cycle loop_nnbrx
2331  END IF
2332  EXIT loop_nnbrx
2333  END DO loop_nnbrx
2334 
2335  istopg = 0
2336 !$OMP PARALLEL DEFAULT(NONE) &
2337 !$OMP PRIVATE(iam,npr,iat,iat1,iat2,lot,istop,tcoord,tcoord2, &
2338 !$OMP tener,tener2,txyz,f2ij,f3ij,f3ik,npjx,npjkx) &
2339 !$OMP SHARED (nat,nnbrx,lsta,lstb,rel,ener,ener2,fxyz,coord,coord2,istopg)
2340 
2341  npr = 1
2342 !$ npr = omp_get_num_threads()
2343  iam = 0
2344 !$ iam = omp_get_thread_num()
2345 
2346  npjx = 300; npjkx = 6000
2347 
2348  IF (npr .NE. 1) THEN
2349 ! PARALLEL CASE
2350 ! create temporary private scalars for reduction sum on energies and
2351 ! temporary private array for reduction sum on forces
2352 !$OMP CRITICAL
2353  ALLOCATE (txyz(3, nat), f2ij(3, npjx), f3ij(3, npjkx), f3ik(3, npjkx))
2354 !$OMP END CRITICAL
2355  IF (iam .EQ. 0) THEN
2356  ener = 0.e0_dp
2357  ener2 = 0.e0_dp
2358  coord = 0.e0_dp
2359  coord2 = 0.e0_dp
2360  END IF
2361 !$OMP DO
2362  DO iat = 1, nat
2363  fxyz(1, iat) = 0.e0_dp
2364  fxyz(2, iat) = 0.e0_dp
2365  fxyz(3, iat) = 0.e0_dp
2366  END DO
2367 !$OMP BARRIER
2368 
2369 ! Each thread treats at most lot atoms
2370  lot = int(real(nat, kind=dp)/real(npr, kind=dp) + .999999999999e0_dp)
2371  iat1 = iam*lot + 1
2372  iat2 = min((iam + 1)*lot, nat)
2373 ! write(*,*) 'subfeniat:iat1,iat2,iam',iat1,iat2,iam
2374  CALL subfeniat_l(iat1, iat2, nat, lsta, lstb, rel, tener, tener2, &
2375  tcoord, tcoord2, nnbrx, txyz, f2ij, npjx, f3ij, npjkx, f3ik, istop)
2376 !$OMP CRITICAL
2377  ener = ener + tener
2378  ener2 = ener2 + tener2
2379  coord = coord + tcoord
2380  coord2 = coord2 + tcoord2
2381  istopg = istopg + istop
2382  DO iat = 1, nat
2383  fxyz(1, iat) = fxyz(1, iat) + txyz(1, iat)
2384  fxyz(2, iat) = fxyz(2, iat) + txyz(2, iat)
2385  fxyz(3, iat) = fxyz(3, iat) + txyz(3, iat)
2386  END DO
2387  DEALLOCATE (txyz, f2ij, f3ij, f3ik)
2388 !$OMP END CRITICAL
2389 
2390  ELSE
2391 ! SERIAL CASE
2392  iat1 = 1
2393  iat2 = nat
2394  ALLOCATE (f2ij(3, npjx), f3ij(3, npjkx), f3ik(3, npjkx))
2395  CALL subfeniat_l(iat1, iat2, nat, lsta, lstb, rel, ener, ener2, &
2396  coord, coord2, nnbrx, fxyz, f2ij, npjx, f3ij, npjkx, f3ik, istopg)
2397  DEALLOCATE (f2ij, f3ij, f3ik)
2398 
2399  END IF
2400 !$OMP END PARALLEL
2401 
2402 ! write(*,*) 'ener,norm force', &
2403 ! ener,DNRM2(3*nat,fxyz,1)
2404  IF (istopg .GT. 0) cpabort("DIMENSION ERROR (see WARNING above)")
2405  ener_var = ener2/nat - (ener/nat)**2
2406  coord = coord/nat
2407  coord_var = coord2/nat - coord**2
2408 
2409  DEALLOCATE (rxyz, icell, lay, lsta, lstb, rel)
2410 
2411  END SUBROUTINE eip_lenosky_silicon
2412 
2413 ! **************************************************************************************************
2414 !> \brief ...
2415 !> \param iat1 ...
2416 !> \param iat2 ...
2417 !> \param nat ...
2418 !> \param lsta ...
2419 !> \param lstb ...
2420 !> \param rel ...
2421 !> \param tener ...
2422 !> \param tener2 ...
2423 !> \param tcoord ...
2424 !> \param tcoord2 ...
2425 !> \param nnbrx ...
2426 !> \param txyz ...
2427 !> \param f2ij ...
2428 !> \param npjx ...
2429 !> \param f3ij ...
2430 !> \param npjkx ...
2431 !> \param f3ik ...
2432 !> \param istop ...
2433 ! **************************************************************************************************
2434  SUBROUTINE subfeniat_l(iat1, iat2, nat, lsta, lstb, rel, tener, tener2, &
2435  tcoord, tcoord2, nnbrx, txyz, f2ij, npjx, f3ij, npjkx, f3ik, istop)
2436 ! for a subset of atoms iat1 to iat2 the routine calculates the (partial) forces
2437 ! txyz acting on these atoms as well as on the atoms (jat, kat) interacting
2438 ! with them and their contribution to the energy (tener).
2439 ! In addition the coordination number tcoord and the second moment of the
2440 ! local energy tener2 and coordination number tcoord2 are returned
2441  INTEGER :: iat1, iat2, nat, lsta, lstb
2442  REAL(kind=dp) :: rel, tener, tener2, tcoord, tcoord2
2443  INTEGER :: nnbrx
2444  REAL(kind=dp) :: txyz, f2ij
2445  INTEGER :: npjx
2446  REAL(kind=dp) :: f3ij
2447  INTEGER :: npjkx
2448  REAL(kind=dp) :: f3ik
2449  INTEGER :: istop
2450 
2451  dimension lsta(2, nat), lstb(nnbrx*nat), rel(5, nnbrx*nat), txyz(3, nat)
2452  dimension f2ij(3, npjx), f3ij(3, npjkx), f3ik(3, npjkx)
2453  REAL(kind=dp), PARAMETER :: tmin_phi = 0.1500000e+01_dp
2454  REAL(kind=dp), PARAMETER :: tmax_phi = 0.4500000e+01_dp
2455  REAL(kind=dp), PARAMETER :: hi_phi = 3.00000000000e0_dp
2456  REAL(kind=dp), PARAMETER :: hsixth_phi = 5.55555555555556e-002_dp
2457  REAL(kind=dp), PARAMETER :: h2sixth_phi = 1.85185185185185e-002_dp
2458  REAL(kind=dp), PARAMETER, DIMENSION(0:9) :: cof_phi = &
2459  (/0.69299400000000e+01_dp, -0.43995000000000e+00_dp, &
2460  -0.17012300000000e+01_dp, -0.16247300000000e+01_dp, &
2461  -0.99696000000000e+00_dp, -0.27391000000000e+00_dp, &
2462  -0.24990000000000e-01_dp, -0.17840000000000e-01_dp, &
2463  -0.96100000000000e-02_dp, 0.00000000000000e+00_dp/)
2464  REAL(kind=dp), PARAMETER, DIMENSION(0:9) :: dof_phi = &
2465  (/0.16533229480429e+03_dp, 0.39415410391417e+02_dp, &
2466  0.68710036300407e+01_dp, 0.53406950884203e+01_dp, &
2467  0.15347960162782e+01_dp, -0.63347591535331e+01_dp, &
2468  -0.17987794021458e+01_dp, 0.47429676211617e+00_dp, &
2469  -0.40087646318907e-01_dp, -0.23942617684055e+00_dp/)
2470  REAL(kind=dp), PARAMETER :: tmin_rho = 0.1500000e+01_dp
2471  REAL(kind=dp), PARAMETER :: tmax_rho = 0.3500000e+01_dp
2472  REAL(kind=dp), PARAMETER :: hi_rho = 5.00000000000e0_dp
2473  REAL(kind=dp), PARAMETER :: hsixth_rho = 3.33333333333333e-002_dp
2474  REAL(kind=dp), PARAMETER :: h2sixth_rho = 6.66666666666667e-003_dp
2475  REAL(kind=dp), PARAMETER, DIMENSION(0:10) :: cof_rho = &
2476  (/0.13747000000000e+00_dp, -0.14831000000000e+00_dp, &
2477  -0.55972000000000e+00_dp, -0.73110000000000e+00_dp, &
2478  -0.76283000000000e+00_dp, -0.72918000000000e+00_dp, &
2479  -0.66620000000000e+00_dp, -0.57328000000000e+00_dp, &
2480  -0.40690000000000e+00_dp, -0.16662000000000e+00_dp, &
2481  0.00000000000000e+00_dp/)
2482  REAL(kind=dp), PARAMETER, DIMENSION(0:10) :: dof_rho = &
2483  (/-0.32275496741918e+01_dp, -0.64119006516165e+01_dp, &
2484  0.10030652280658e+02_dp, 0.22937915289857e+01_dp, &
2485  0.17416816033995e+01_dp, 0.54648205741626e+00_dp, &
2486  0.47189016693543e+00_dp, 0.20569572748420e+01_dp, &
2487  0.23192807336964e+01_dp, -0.24908020962757e+00_dp, &
2488  -0.12371959895186e+02_dp/)
2489  REAL(kind=dp), PARAMETER :: tmin_fff = 0.1500000e+01_dp
2490  REAL(kind=dp), PARAMETER :: tmax_fff = 0.3500000e+01_dp
2491  REAL(kind=dp), PARAMETER :: hi_fff = 4.50000000000e0_dp
2492  REAL(kind=dp), PARAMETER :: hsixth_fff = 3.70370370370370e-002_dp
2493  REAL(kind=dp), PARAMETER :: h2sixth_fff = 8.23045267489712e-003_dp
2494  REAL(kind=dp), PARAMETER, DIMENSION(0:9) :: cof_fff = &
2495  (/0.12503100000000e+01_dp, 0.86821000000000e+00_dp, &
2496  0.60846000000000e+00_dp, 0.48756000000000e+00_dp, &
2497  0.44163000000000e+00_dp, 0.37610000000000e+00_dp, &
2498  0.27145000000000e+00_dp, 0.14814000000000e+00_dp, &
2499  0.48550000000000e-01_dp, 0.00000000000000e+00_dp/)
2500  REAL(kind=dp), PARAMETER, DIMENSION(0:9) :: dof_fff = &
2501  (/0.27904652711432e+02_dp, -0.45230754228635e+01_dp, &
2502  0.50531739800222e+01_dp, 0.11806545027747e+01_dp, &
2503  -0.66693699112098e+00_dp, -0.89430653829079e+00_dp, &
2504  -0.50891685571587e+00_dp, 0.66278396115427e+00_dp, &
2505  0.73976101109878e+00_dp, 0.25795319944506e+01_dp/)
2506  REAL(kind=dp), PARAMETER :: tmin_uuu = -0.1770930e+01_dp
2507  REAL(kind=dp), PARAMETER :: tmax_uuu = 0.7908520e+01_dp
2508  REAL(kind=dp), PARAMETER :: hi_uuu = 0.723181585730594e0_dp
2509  REAL(kind=dp), PARAMETER :: hsixth_uuu = 0.230463095238095e0_dp
2510  REAL(kind=dp), PARAMETER :: h2sixth_uuu = 0.318679429600340e0_dp
2511  REAL(kind=dp), PARAMETER, DIMENSION(0:7) :: cof_uuu = &
2512  (/-0.10749300000000e+01_dp, -0.20045000000000e+00_dp, &
2513  0.41422000000000e+00_dp, 0.87939000000000e+00_dp, &
2514  0.12668900000000e+01_dp, 0.16299800000000e+01_dp, &
2515  0.19773800000000e+01_dp, 0.23961800000000e+01_dp/)
2516  REAL(kind=dp), PARAMETER, DIMENSION(0:7) :: dof_uuu = &
2517  (/-0.14827125747284e+00_dp, -0.14922155328475e+00_dp, &
2518  -0.70113224223509e-01_dp, -0.39449020349230e-01_dp, &
2519  -0.15815242579643e-01_dp, 0.26112640061855e-01_dp, &
2520  -0.13786974745095e+00_dp, 0.74941595372657e+00_dp/)
2521  REAL(kind=dp), PARAMETER :: tmin_ggg = -0.1000000e+01_dp
2522  REAL(kind=dp), PARAMETER :: tmax_ggg = 0.8001400e+00_dp
2523  REAL(kind=dp), PARAMETER :: hi_ggg = 3.88858644327663e0_dp
2524  REAL(kind=dp), PARAMETER :: hsixth_ggg = 4.28604761904762e-002_dp
2525  REAL(kind=dp), PARAMETER :: h2sixth_ggg = 1.10221225156463e-002_dp
2526  REAL(kind=dp), PARAMETER, DIMENSION(0:7) :: cof_ggg = &
2527  (/0.52541600000000e+01_dp, 0.23591500000000e+01_dp, &
2528  0.11959500000000e+01_dp, 0.12299500000000e+01_dp, &
2529  0.20356500000000e+01_dp, 0.34247400000000e+01_dp, &
2530  0.49485900000000e+01_dp, 0.56179900000000e+01_dp/)
2531  REAL(kind=dp), PARAMETER, DIMENSION(0:7) :: dof_ggg = &
2532  (/0.15826876132396e+02_dp, 0.31176239377907e+02_dp, &
2533  0.16589446539683e+02_dp, 0.11083892500520e+02_dp, &
2534  0.90887216383860e+01_dp, 0.54902279653967e+01_dp, &
2535  -0.18823313223755e+02_dp, -0.77183416481005e+01_dp/)
2536 
2537  REAL(kind=dp) :: a2_fff, a2_ggg, a_fff, a_ggg, b2_fff, b2_ggg, b_fff, &
2538  b_ggg, cof1_fff, cof1_ggg, cof2_fff, cof2_ggg, cof3_fff, &
2539  cof3_ggg, cof4_fff, cof4_ggg, cof_fff_khi, cof_fff_klo, &
2540  cof_ggg_khi, cof_ggg_klo, coord_iat, costheta, dens, &
2541  dens2, dens3, dof_fff_khi, dof_fff_klo, dof_ggg_khi, &
2542  dof_ggg_klo, e_phi, e_uuu, ener_iat, ep_phi, ep_uuu, &
2543  fij, fijp, fik, fikp, fxij, fxik, fyij, fyik, fzij, fzik, &
2544  gjik, gjikp, rho, rhop, rij, rik, sij, sik, t1, t2, t3, t4, &
2545  tt, tt_fff, tt_ggg, xarg, ypt1_fff, ypt1_ggg, ypt2_fff, &
2546  ypt2_ggg, yt1_fff, yt1_ggg, yt2_fff, yt2_ggg
2547 
2548  INTEGER :: iat, jat, jbr, jcnt, jkcnt, kat, kbr, khi_fff, khi_ggg, &
2549  klo_fff, klo_ggg
2550 
2551 ! initialize temporary private scalars for reduction sum on energies and
2552 ! private workarray txyz for forces forces
2553  tener = 0.e0_dp
2554  tener2 = 0.e0_dp
2555  tcoord = 0.e0_dp
2556  tcoord2 = 0.e0_dp
2557  istop = 0
2558  DO iat = 1, nat
2559  txyz(1, iat) = 0.e0_dp
2560  txyz(2, iat) = 0.e0_dp
2561  txyz(3, iat) = 0.e0_dp
2562  END DO
2563 
2564 ! calculation of forces, energy
2565 
2566  forces_and_energy: DO iat = iat1, iat2
2567 
2568  dens2 = 0.e0_dp
2569  dens3 = 0.e0_dp
2570  jcnt = 0
2571  jkcnt = 0
2572  coord_iat = 0.e0_dp
2573  ener_iat = 0.e0_dp
2574  calculate: DO jbr = lsta(1, iat), lsta(2, iat)
2575  jat = lstb(jbr)
2576  jcnt = jcnt + 1
2577  IF (jcnt .GT. npjx) THEN
2578  WRITE (*, *) 'WARNING: enlarge npjx'
2579  istop = 1
2580  RETURN
2581  END IF
2582 
2583  fxij = rel(1, jbr)
2584  fyij = rel(2, jbr)
2585  fzij = rel(3, jbr)
2586  rij = rel(4, jbr)
2587  sij = rel(5, jbr)
2588 
2589 ! coordination number calculated with soft cutoff between first
2590 ! nearest neighbor and midpoint of first and second nearest neighbor
2591  IF (rij .LE. 2.36e0_dp) THEN
2592  coord_iat = coord_iat + 1.e0_dp
2593  ELSE IF (rij .GE. 3.12e0_dp) THEN
2594  ELSE
2595  xarg = (rij - 2.36e0_dp)*(1.e0_dp/(3.12e0_dp - 2.36e0_dp))
2596  coord_iat = coord_iat + (2*xarg + 1.e0_dp)*(xarg - 1.e0_dp)**2
2597  END IF
2598 
2599 ! pairpotential term
2600  CALL splint(cof_phi, dof_phi, tmin_phi, tmax_phi, &
2601  hsixth_phi, h2sixth_phi, hi_phi, 10, rij, e_phi, ep_phi)
2602  ener_iat = ener_iat + (e_phi*.5e0_dp)
2603  txyz(1, iat) = txyz(1, iat) - fxij*(ep_phi*.5e0_dp)
2604  txyz(2, iat) = txyz(2, iat) - fyij*(ep_phi*.5e0_dp)
2605  txyz(3, iat) = txyz(3, iat) - fzij*(ep_phi*.5e0_dp)
2606  txyz(1, jat) = txyz(1, jat) + fxij*(ep_phi*.5e0_dp)
2607  txyz(2, jat) = txyz(2, jat) + fyij*(ep_phi*.5e0_dp)
2608  txyz(3, jat) = txyz(3, jat) + fzij*(ep_phi*.5e0_dp)
2609 
2610 ! 2 body embedding term
2611  CALL splint(cof_rho, dof_rho, tmin_rho, tmax_rho, &
2612  hsixth_rho, h2sixth_rho, hi_rho, 11, rij, rho, rhop)
2613  dens2 = dens2 + rho
2614  f2ij(1, jcnt) = fxij*rhop
2615  f2ij(2, jcnt) = fyij*rhop
2616  f2ij(3, jcnt) = fzij*rhop
2617 
2618 ! 3 body embedding term
2619  CALL splint(cof_fff, dof_fff, tmin_fff, tmax_fff, &
2620  hsixth_fff, h2sixth_fff, hi_fff, 10, rij, fij, fijp)
2621 
2622  embed_3body: DO kbr = lsta(1, iat), lsta(2, iat)
2623  kat = lstb(kbr)
2624  IF (kat .LT. jat) THEN
2625  jkcnt = jkcnt + 1
2626  IF (jkcnt .GT. npjkx) THEN
2627  WRITE (*, *) 'WARNING: enlarge npjkx', npjkx
2628  istop = 1
2629  RETURN
2630  END IF
2631 
2632 ! begin unoptimized original version:
2633 ! fxik=rel(1,kbr)
2634 ! fyik=rel(2,kbr)
2635 ! fzik=rel(3,kbr)
2636 ! rik=rel(4,kbr)
2637 ! sik=rel(5,kbr)
2638 !
2639 ! call splint(cof_fff,dof_fff,tmin_fff,tmax_fff, &
2640 ! hsixth_fff,h2sixth_fff,hi_fff,10,rik,fik,fikp)
2641 ! costheta=fxij*fxik+fyij*fyik+fzij*fzik
2642 ! call splint(cof_ggg,dof_ggg,tmin_ggg,tmax_ggg, &
2643 ! hsixth_ggg,h2sixth_ggg,hi_ggg,8,costheta,gjik,gjikp)
2644 ! end unoptimized original version:
2645 
2646 ! begin optimized version
2647  rik = rel(4, kbr)
2648  IF (rik .GT. tmax_fff) THEN
2649  fikp = 0.e0_dp; fik = 0.e0_dp
2650  gjik = 0.e0_dp; gjikp = 0.e0_dp; sik = 0.e0_dp
2651  costheta = 0.e0_dp; fxik = 0.e0_dp; fyik = 0.e0_dp; fzik = 0.e0_dp
2652  ELSE IF (rik .LT. tmin_fff) THEN
2653  fxik = rel(1, kbr)
2654  fyik = rel(2, kbr)
2655  fzik = rel(3, kbr)
2656  costheta = fxij*fxik + fyij*fyik + fzij*fzik
2657  sik = rel(5, kbr)
2658  fikp = hi_fff*(cof_fff(1) - cof_fff(0)) - &
2659  (dof_fff(1) + 2.e0_dp*dof_fff(0))*hsixth_fff
2660  fik = cof_fff(0) + (rik - tmin_fff)*fikp
2661  tt_ggg = (costheta - tmin_ggg)*hi_ggg
2662  IF (costheta .GT. tmax_ggg) THEN
2663  gjikp = hi_ggg*(cof_ggg(8 - 1) - cof_ggg(8 - 2)) + &
2664  (2.e0_dp*dof_ggg(8 - 1) + dof_ggg(8 - 2))*hsixth_ggg
2665  gjik = cof_ggg(8 - 1) + (costheta - tmax_ggg)*gjikp
2666  ELSE
2667  klo_ggg = int(tt_ggg)
2668  khi_ggg = klo_ggg + 1
2669  cof_ggg_klo = cof_ggg(klo_ggg)
2670  dof_ggg_klo = dof_ggg(klo_ggg)
2671  b_ggg = tt_ggg - klo_ggg
2672  a_ggg = 1.e0_dp - b_ggg
2673  cof_ggg_khi = cof_ggg(khi_ggg)
2674  dof_ggg_khi = dof_ggg(khi_ggg)
2675  b2_ggg = b_ggg*b_ggg
2676  gjik = a_ggg*cof_ggg_klo
2677  gjikp = cof_ggg_khi - cof_ggg_klo
2678  a2_ggg = a_ggg*a_ggg
2679  cof1_ggg = a2_ggg - 1.e0_dp
2680  cof2_ggg = b2_ggg - 1.e0_dp
2681  gjik = gjik + b_ggg*cof_ggg_khi
2682  gjikp = hi_ggg*gjikp
2683  cof3_ggg = 3.e0_dp*b2_ggg
2684  cof4_ggg = 3.e0_dp*a2_ggg
2685  cof1_ggg = a_ggg*cof1_ggg
2686  cof2_ggg = b_ggg*cof2_ggg
2687  cof3_ggg = cof3_ggg - 1.e0_dp
2688  cof4_ggg = cof4_ggg - 1.e0_dp
2689  yt1_ggg = cof1_ggg*dof_ggg_klo
2690  yt2_ggg = cof2_ggg*dof_ggg_khi
2691  ypt1_ggg = cof3_ggg*dof_ggg_khi
2692  ypt2_ggg = cof4_ggg*dof_ggg_klo
2693  gjik = gjik + (yt1_ggg + yt2_ggg)*h2sixth_ggg
2694  gjikp = gjikp + (ypt1_ggg - ypt2_ggg)*hsixth_ggg
2695  END IF
2696  ELSE
2697  fxik = rel(1, kbr)
2698  tt_fff = rik - tmin_fff
2699  costheta = fxij*fxik
2700  fyik = rel(2, kbr)
2701  tt_fff = tt_fff*hi_fff
2702  costheta = costheta + fyij*fyik
2703  fzik = rel(3, kbr)
2704  klo_fff = int(tt_fff)
2705  costheta = costheta + fzij*fzik
2706  sik = rel(5, kbr)
2707  tt_ggg = (costheta - tmin_ggg)*hi_ggg
2708  IF (costheta .GT. tmax_ggg) THEN
2709  gjikp = hi_ggg*(cof_ggg(8 - 1) - cof_ggg(8 - 2)) + &
2710  (2.e0_dp*dof_ggg(8 - 1) + dof_ggg(8 - 2))*hsixth_ggg
2711  gjik = cof_ggg(8 - 1) + (costheta - tmax_ggg)*gjikp
2712  khi_fff = klo_fff + 1
2713  cof_fff_klo = cof_fff(klo_fff)
2714  dof_fff_klo = dof_fff(klo_fff)
2715  b_fff = tt_fff - klo_fff
2716  a_fff = 1.e0_dp - b_fff
2717  cof_fff_khi = cof_fff(khi_fff)
2718  dof_fff_khi = dof_fff(khi_fff)
2719  b2_fff = b_fff*b_fff
2720  fik = a_fff*cof_fff_klo
2721  fikp = cof_fff_khi - cof_fff_klo
2722  a2_fff = a_fff*a_fff
2723  cof1_fff = a2_fff - 1.e0_dp
2724  cof2_fff = b2_fff - 1.e0_dp
2725  fik = fik + b_fff*cof_fff_khi
2726  fikp = hi_fff*fikp
2727  cof3_fff = 3.e0_dp*b2_fff
2728  cof4_fff = 3.e0_dp*a2_fff
2729  cof1_fff = a_fff*cof1_fff
2730  cof2_fff = b_fff*cof2_fff
2731  cof3_fff = cof3_fff - 1.e0_dp
2732  cof4_fff = cof4_fff - 1.e0_dp
2733  yt1_fff = cof1_fff*dof_fff_klo
2734  yt2_fff = cof2_fff*dof_fff_khi
2735  ypt1_fff = cof3_fff*dof_fff_khi
2736  ypt2_fff = cof4_fff*dof_fff_klo
2737  fik = fik + (yt1_fff + yt2_fff)*h2sixth_fff
2738  fikp = fikp + (ypt1_fff - ypt2_fff)*hsixth_fff
2739  ELSE
2740  klo_ggg = int(tt_ggg)
2741  khi_ggg = klo_ggg + 1
2742  khi_fff = klo_fff + 1
2743  cof_ggg_klo = cof_ggg(klo_ggg)
2744  cof_fff_klo = cof_fff(klo_fff)
2745  dof_ggg_klo = dof_ggg(klo_ggg)
2746  dof_fff_klo = dof_fff(klo_fff)
2747  b_ggg = tt_ggg - klo_ggg
2748  b_fff = tt_fff - klo_fff
2749  a_ggg = 1.e0_dp - b_ggg
2750  a_fff = 1.e0_dp - b_fff
2751  cof_ggg_khi = cof_ggg(khi_ggg)
2752  cof_fff_khi = cof_fff(khi_fff)
2753  dof_ggg_khi = dof_ggg(khi_ggg)
2754  dof_fff_khi = dof_fff(khi_fff)
2755  b2_ggg = b_ggg*b_ggg
2756  b2_fff = b_fff*b_fff
2757  gjik = a_ggg*cof_ggg_klo
2758  fik = a_fff*cof_fff_klo
2759  gjikp = cof_ggg_khi - cof_ggg_klo
2760  fikp = cof_fff_khi - cof_fff_klo
2761  a2_ggg = a_ggg*a_ggg
2762  a2_fff = a_fff*a_fff
2763  cof1_ggg = a2_ggg - 1.e0_dp
2764  cof1_fff = a2_fff - 1.e0_dp
2765  cof2_ggg = b2_ggg - 1.e0_dp
2766  cof2_fff = b2_fff - 1.e0_dp
2767  gjik = gjik + b_ggg*cof_ggg_khi
2768  fik = fik + b_fff*cof_fff_khi
2769  gjikp = hi_ggg*gjikp
2770  fikp = hi_fff*fikp
2771  cof3_ggg = 3.e0_dp*b2_ggg
2772  cof3_fff = 3.e0_dp*b2_fff
2773  cof4_ggg = 3.e0_dp*a2_ggg
2774  cof4_fff = 3.e0_dp*a2_fff
2775  cof1_ggg = a_ggg*cof1_ggg
2776  cof1_fff = a_fff*cof1_fff
2777  cof2_ggg = b_ggg*cof2_ggg
2778  cof2_fff = b_fff*cof2_fff
2779  cof3_ggg = cof3_ggg - 1.e0_dp
2780  cof3_fff = cof3_fff - 1.e0_dp
2781  cof4_ggg = cof4_ggg - 1.e0_dp
2782  cof4_fff = cof4_fff - 1.e0_dp
2783  yt1_ggg = cof1_ggg*dof_ggg_klo
2784  yt1_fff = cof1_fff*dof_fff_klo
2785  yt2_ggg = cof2_ggg*dof_ggg_khi
2786  yt2_fff = cof2_fff*dof_fff_khi
2787  ypt1_ggg = cof3_ggg*dof_ggg_khi
2788  ypt1_fff = cof3_fff*dof_fff_khi
2789  ypt2_ggg = cof4_ggg*dof_ggg_klo
2790  ypt2_fff = cof4_fff*dof_fff_klo
2791  gjik = gjik + (yt1_ggg + yt2_ggg)*h2sixth_ggg
2792  fik = fik + (yt1_fff + yt2_fff)*h2sixth_fff
2793  gjikp = gjikp + (ypt1_ggg - ypt2_ggg)*hsixth_ggg
2794  fikp = fikp + (ypt1_fff - ypt2_fff)*hsixth_fff
2795  END IF
2796  END IF
2797 ! end optimized version
2798 
2799  tt = fij*fik
2800  dens3 = dens3 + tt*gjik
2801 
2802  t1 = fijp*fik*gjik
2803  t2 = sij*(tt*gjikp)
2804  f3ij(1, jkcnt) = fxij*t1 + (fxik - fxij*costheta)*t2
2805  f3ij(2, jkcnt) = fyij*t1 + (fyik - fyij*costheta)*t2
2806  f3ij(3, jkcnt) = fzij*t1 + (fzik - fzij*costheta)*t2
2807 
2808  t3 = fikp*fij*gjik
2809  t4 = sik*(tt*gjikp)
2810  f3ik(1, jkcnt) = fxik*t3 + (fxij - fxik*costheta)*t4
2811  f3ik(2, jkcnt) = fyik*t3 + (fyij - fyik*costheta)*t4
2812  f3ik(3, jkcnt) = fzik*t3 + (fzij - fzik*costheta)*t4
2813  END IF
2814 
2815  END DO embed_3body
2816  END DO calculate
2817 
2818  dens = dens2 + dens3
2819  CALL splint(cof_uuu, dof_uuu, tmin_uuu, tmax_uuu, &
2820  hsixth_uuu, h2sixth_uuu, hi_uuu, 8, dens, e_uuu, ep_uuu)
2821  ener_iat = ener_iat + e_uuu
2822 
2823 ! Only now ep_uu is known and the forces can be calculated, lets loop again
2824  jcnt = 0
2825  jkcnt = 0
2826  loop_again: DO jbr = lsta(1, iat), lsta(2, iat)
2827  jat = lstb(jbr)
2828  jcnt = jcnt + 1
2829  txyz(1, iat) = txyz(1, iat) - ep_uuu*f2ij(1, jcnt)
2830  txyz(2, iat) = txyz(2, iat) - ep_uuu*f2ij(2, jcnt)
2831  txyz(3, iat) = txyz(3, iat) - ep_uuu*f2ij(3, jcnt)
2832  txyz(1, jat) = txyz(1, jat) + ep_uuu*f2ij(1, jcnt)
2833  txyz(2, jat) = txyz(2, jat) + ep_uuu*f2ij(2, jcnt)
2834  txyz(3, jat) = txyz(3, jat) + ep_uuu*f2ij(3, jcnt)
2835 
2836 ! 3 body embedding term
2837  DO kbr = lsta(1, iat), lsta(2, iat)
2838  kat = lstb(kbr)
2839  IF (kat .LT. jat) THEN
2840  jkcnt = jkcnt + 1
2841 
2842  txyz(1, iat) = txyz(1, iat) - ep_uuu*(f3ij(1, jkcnt) + f3ik(1, jkcnt))
2843  txyz(2, iat) = txyz(2, iat) - ep_uuu*(f3ij(2, jkcnt) + f3ik(2, jkcnt))
2844  txyz(3, iat) = txyz(3, iat) - ep_uuu*(f3ij(3, jkcnt) + f3ik(3, jkcnt))
2845  txyz(1, jat) = txyz(1, jat) + ep_uuu*f3ij(1, jkcnt)
2846  txyz(2, jat) = txyz(2, jat) + ep_uuu*f3ij(2, jkcnt)
2847  txyz(3, jat) = txyz(3, jat) + ep_uuu*f3ij(3, jkcnt)
2848  txyz(1, kat) = txyz(1, kat) + ep_uuu*f3ik(1, jkcnt)
2849  txyz(2, kat) = txyz(2, kat) + ep_uuu*f3ik(2, jkcnt)
2850  txyz(3, kat) = txyz(3, kat) + ep_uuu*f3ik(3, jkcnt)
2851  END IF
2852  END DO
2853 
2854  END DO loop_again
2855 
2856 ! write(*,'(a,i4,x,e19.12,x,e10.3)') 'iat,ener_iat,coord_iat', &
2857 ! iat,ener_iat,coord_iat
2858  tener = tener + ener_iat
2859  tener2 = tener2 + ener_iat**2
2860  tcoord = tcoord + coord_iat
2861  tcoord2 = tcoord2 + coord_iat**2
2862 
2863  END DO forces_and_energy
2864 
2865  END SUBROUTINE subfeniat_l
2866 
2867 ! **************************************************************************************************
2868 !> \brief ...
2869 !> \param iat ...
2870 !> \param nn ...
2871 !> \param ncx ...
2872 !> \param ll1 ...
2873 !> \param ll2 ...
2874 !> \param ll3 ...
2875 !> \param l1 ...
2876 !> \param l2 ...
2877 !> \param l3 ...
2878 !> \param myspace ...
2879 !> \param rxyz ...
2880 !> \param icell ...
2881 !> \param lstb ...
2882 !> \param lay ...
2883 !> \param rel ...
2884 !> \param cut2 ...
2885 !> \param indlst ...
2886 ! **************************************************************************************************
2887  SUBROUTINE sublstiat_l(iat, nn, ncx, ll1, ll2, ll3, l1, l2, l3, myspace, &
2888  rxyz, icell, lstb, lay, rel, cut2, indlst)
2889 ! finds the neighbours of atom iat (specified by lsta and lstb) and and
2890 ! the relative position rel of iat with respect to these neighbours
2891  INTEGER :: iat, nn, ncx, ll1, ll2, ll3, l1, l2, l3, &
2892  myspace
2893  REAL(kind=dp) :: rxyz
2894  INTEGER :: icell, lstb, lay
2895  REAL(kind=dp) :: rel, cut2
2896  INTEGER :: indlst
2897 
2898  dimension rxyz(3, nn), lay(nn), icell(0:ncx, -1:ll1, -1:ll2, -1:ll3), &
2899  lstb(0:myspace - 1), rel(5, 0:myspace - 1)
2900 
2901  INTEGER :: jat, jj, k1, k2, k3
2902  REAL(kind=dp) :: rr2, tt, xrel, yrel, zrel, tti
2903 
2904  loop_k3: DO k3 = l3 - 1, l3 + 1
2905  loop_k2: DO k2 = l2 - 1, l2 + 1
2906  loop_k1: DO k1 = l1 - 1, l1 + 1
2907  loop_jj: DO jj = 1, icell(0, k1, k2, k3)
2908  jat = icell(jj, k1, k2, k3)
2909  IF (jat .EQ. iat) cycle loop_k3
2910  xrel = rxyz(1, iat) - rxyz(1, jat)
2911  yrel = rxyz(2, iat) - rxyz(2, jat)
2912  zrel = rxyz(3, iat) - rxyz(3, jat)
2913  rr2 = xrel**2 + yrel**2 + zrel**2
2914  IF (rr2 .LE. cut2) THEN
2915  indlst = min(indlst, myspace - 1)
2916  lstb(indlst) = lay(jat)
2917 ! write(*,*) 'iat,indlst,lay(jat)',iat,indlst,lay(jat)
2918  tt = sqrt(rr2)
2919  tti = 1.e0_dp/tt
2920  rel(1, indlst) = xrel*tti
2921  rel(2, indlst) = yrel*tti
2922  rel(3, indlst) = zrel*tti
2923  rel(4, indlst) = tt
2924  rel(5, indlst) = tti
2925  indlst = indlst + 1
2926  END IF
2927  END DO loop_jj
2928  END DO loop_k1
2929  END DO loop_k2
2930  END DO loop_k3
2931 
2932  RETURN
2933  END SUBROUTINE sublstiat_l
2934 
2935 ! **************************************************************************************************
2936 !> \brief ...
2937 !> \param ya ...
2938 !> \param y2a ...
2939 !> \param tmin ...
2940 !> \param tmax ...
2941 !> \param hsixth ...
2942 !> \param h2sixth ...
2943 !> \param hi ...
2944 !> \param n ...
2945 !> \param x ...
2946 !> \param y ...
2947 !> \param yp ...
2948 ! **************************************************************************************************
2949  SUBROUTINE splint(ya, y2a, tmin, tmax, hsixth, h2sixth, hi, n, x, y, yp)
2950  REAL(kind=dp) :: ya, y2a, tmin, tmax, hsixth, h2sixth, hi
2951  INTEGER :: n
2952  REAL(kind=dp) :: x, y, yp
2953 
2954  dimension y2a(0:n - 1), ya(0:n - 1)
2955  REAL(kind=dp) :: a, a2, b, b2, cof1, cof2, cof3, cof4, tt, &
2956  y2a_khi, ya_klo, y2a_klo, ya_khi, ypt1, ypt2, yt1, yt2
2957  INTEGER :: klo, khi
2958 
2959 ! interpolate if the argument is outside the cubic spline interval [tmin,tmax]
2960  tt = (x - tmin)*hi
2961  IF (x .LT. tmin) THEN
2962  yp = hi*(ya(1) - ya(0)) - &
2963  (y2a(1) + 2.e0_dp*y2a(0))*hsixth
2964  y = ya(0) + (x - tmin)*yp
2965  ELSE IF (x .GT. tmax) THEN
2966  yp = hi*(ya(n - 1) - ya(n - 2)) + &
2967  (2.e0_dp*y2a(n - 1) + y2a(n - 2))*hsixth
2968  y = ya(n - 1) + (x - tmax)*yp
2969 ! otherwise evaluate cubic spline
2970  ELSE
2971  klo = int(tt)
2972  khi = klo + 1
2973  ya_klo = ya(klo)
2974  y2a_klo = y2a(klo)
2975  b = tt - klo
2976  a = 1.e0_dp - b
2977  ya_khi = ya(khi)
2978  y2a_khi = y2a(khi)
2979  b2 = b*b
2980  y = a*ya_klo
2981  yp = ya_khi - ya_klo
2982  a2 = a*a
2983  cof1 = a2 - 1.e0_dp
2984  cof2 = b2 - 1.e0_dp
2985  y = y + b*ya_khi
2986  yp = hi*yp
2987  cof3 = 3.e0_dp*b2
2988  cof4 = 3.e0_dp*a2
2989  cof1 = a*cof1
2990  cof2 = b*cof2
2991  cof3 = cof3 - 1.e0_dp
2992  cof4 = cof4 - 1.e0_dp
2993  yt1 = cof1*y2a_klo
2994  yt2 = cof2*y2a_khi
2995  ypt1 = cof3*y2a_khi
2996  ypt2 = cof4*y2a_klo
2997  y = y + (yt1 + yt2)*h2sixth
2998  yp = yp + (ypt1 - ypt2)*hsixth
2999  END IF
3000  RETURN
3001  END SUBROUTINE splint
3002 
3003 END 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:117
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:195
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:803
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:1013
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)
returns information about various attributes of the given subsys
Definition: cp_subsys_types.F:345
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_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_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:724
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