(git:374b731)
Loading...
Searching...
No Matches
eip_silicon.F
Go to the documentation of this file.
1!--------------------------------------------------------------------------------------------------!
2! CP2K: A general program to perform molecular dynamics simulations !
3! Copyright 2000-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! **************************************************************************************************
17MODULE eip_silicon
18 USE atomic_kind_list_types, ONLY: atomic_kind_list_type
19 USE atomic_kind_types, ONLY: atomic_kind_type,&
20 get_atomic_kind
21 USE cell_types, ONLY: cell_type,&
22 get_cell
23 USE cp_log_handling, ONLY: cp_get_default_logger,&
24 cp_logger_type
25 USE cp_output_handling, ONLY: cp_p_file,&
26 cp_print_key_finished_output,&
27 cp_print_key_should_output,&
28 cp_print_key_unit_nr
29 USE cp_subsys_types, ONLY: cp_subsys_get,&
30 cp_subsys_type
31 USE distribution_1d_types, ONLY: distribution_1d_type
32 USE eip_environment_types, ONLY: eip_env_get,&
33 eip_environment_type
34 USE input_section_types, ONLY: section_vals_get_subs_vals,&
35 section_vals_type
36 USE kinds, ONLY: dp
37 USE message_passing, ONLY: mp_para_env_type
38 USE particle_types, ONLY: particle_type
39 USE physcon, ONLY: angstrom,&
40 evolt
41
42!$ USE OMP_LIB, ONLY: omp_get_max_threads, omp_get_thread_num, omp_get_num_threads
43#include "./base/base_uses.f90"
44
45 IMPLICIT NONE
46 PRIVATE
47
48 CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'eip_silicon'
49
50 ! *** Public subroutines ***
51 PUBLIC :: eip_bazant, eip_lenosky
52
53!***
54
55CONTAINS
56
57! **************************************************************************************************
58!> \brief Interface routine of Goedecker's Bazant EDIP to CP2K
59!> \param eip_env ...
60!> \par Literature
61!> http://www-math.mit.edu/~bazant/EDIP
62!> M.Z. Bazant & E. Kaxiras: Modeling of Covalent Bonding in Solids by
63!> Inversion of Cohesive Energy Curves;
64!> Phys. Rev. Lett. 77, 4370 (1996)
65!> M.Z. Bazant, E. Kaxiras and J.F. Justo: Environment-dependent interatomic
66!> potential for bulk silicon;
67!> Phys. Rev. B 56, 8542-8552 (1997)
68!> S. Goedecker: Optimization and parallelization of a force field for silicon
69!> using OpenMP; CPC 148, 1 (2002)
70!> \par History
71!> 03.2006 initial create [tdk]
72!> \author Thomas D. Kuehne (tkuehne@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
3003END 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
atomic_kind_list_types::atomic_kind_list_type
represent a list of objects
Definition atomic_kind_list_types.F:45
atomic_kind_types::atomic_kind_type
Provides all information about an atomic kind.
Definition atomic_kind_types.F:45
cell_types::cell_type
Type defining parameters related to the simulation cell.
Definition cell_types.F:55
cp_log_handling::cp_logger_type
type of a logger, at the moment it contains just a print level starting at which level it should be l...
Definition cp_log_handling.F:140
cp_subsys_types::cp_subsys_type
represents a system: atoms, molecules, their pos,vel,...
Definition cp_subsys_types.F:89
distribution_1d_types::distribution_1d_type
structure to store local (to a processor) ordered lists of integers.
Definition distribution_1d_types.F:62
eip_environment_types::eip_environment_type
The empirical interatomic potential environment.
Definition eip_environment_types.F:83
input_section_types::section_vals_type
stores the values of a section
Definition input_section_types.F:126
message_passing::mp_para_env_type
stores all the informations relevant to an mpi environment
Definition message_passing.F:763
particle_types::particle_type
Definition particle_types.F:35