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