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