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