35 #include "./base/base_uses.f90"
41 CHARACTER(len=*),
PARAMETER,
PRIVATE :: moduleN =
'force_env_utils'
67 SUBROUTINE force_env_shake(force_env, dt, shake_tol, log_unit, lagrange_mult, dump_lm, &
68 pos, vel, compold, reset)
70 TYPE(force_env_type),
POINTER :: force_env
71 REAL(kind=
dp),
INTENT(IN),
OPTIONAL :: dt
72 REAL(kind=
dp),
INTENT(IN) :: shake_tol
73 INTEGER,
INTENT(in),
OPTIONAL :: log_unit, lagrange_mult
74 LOGICAL,
INTENT(IN),
OPTIONAL :: dump_lm
75 REAL(kind=
dp),
DIMENSION(:, :),
INTENT(INOUT), &
76 OPTIONAL,
TARGET :: pos, vel
77 LOGICAL,
INTENT(IN),
OPTIONAL :: compold, reset
79 CHARACTER(len=*),
PARAMETER :: routinen =
'force_env_shake'
81 INTEGER :: handle, i, iparticle, iparticle_kind, iparticle_local, j, my_lagrange_mult, &
82 my_log_unit, nparticle_kind, nparticle_local
83 LOGICAL :: has_pos, has_vel, my_dump_lm
85 REAL(kind=
dp),
DIMENSION(:, :),
POINTER :: my_pos, my_vel
86 TYPE(atomic_kind_list_type),
POINTER :: atomic_kinds
87 TYPE(cell_type),
POINTER :: cell
88 TYPE(cp_subsys_type),
POINTER :: subsys
89 TYPE(distribution_1d_type),
POINTER :: local_molecules, local_particles
90 TYPE(global_constraint_type),
POINTER :: gci
91 TYPE(molecule_kind_list_type),
POINTER :: molecule_kinds
92 TYPE(molecule_list_type),
POINTER :: molecules
93 TYPE(particle_list_type),
POINTER :: particles
95 CALL timeset(routinen, handle)
96 cpassert(
ASSOCIATED(force_env))
97 cpassert(force_env%ref_count > 0)
99 IF (
PRESENT(log_unit)) my_log_unit = log_unit
100 my_lagrange_mult = -1
101 IF (
PRESENT(lagrange_mult)) my_lagrange_mult = lagrange_mult
103 IF (
PRESENT(dump_lm)) my_dump_lm = dump_lm
104 NULLIFY (subsys, cell, molecules, molecule_kinds, local_molecules, particles, &
106 IF (
PRESENT(pos)) my_pos => pos
107 IF (
PRESENT(vel)) my_vel => vel
109 IF (
PRESENT(dt)) mydt = dt
112 atomic_kinds=atomic_kinds, &
113 local_molecules=local_molecules, &
114 local_particles=local_particles, &
115 molecules=molecules, &
116 molecule_kinds=molecule_kinds, &
117 particles=particles, &
119 nparticle_kind = atomic_kinds%n_els
120 IF (
PRESENT(compold))
THEN
122 CALL getold(gci, local_molecules, molecules%els, molecule_kinds%els, &
127 IF (.NOT.
ASSOCIATED(my_pos))
THEN
129 ALLOCATE (my_pos(3, particles%n_els))
131 DO iparticle_kind = 1, nparticle_kind
132 nparticle_local = local_particles%n_el(iparticle_kind)
133 DO iparticle_local = 1, nparticle_local
134 iparticle = local_particles%list(iparticle_kind)%array(iparticle_local)
135 my_pos(:, iparticle) = particles%els(iparticle)%r(:)
140 IF (.NOT.
ASSOCIATED(my_vel))
THEN
142 ALLOCATE (my_vel(3, particles%n_els))
144 DO iparticle_kind = 1, nparticle_kind
145 nparticle_local = local_particles%n_el(iparticle_kind)
146 DO iparticle_local = 1, nparticle_local
147 iparticle = local_particles%list(iparticle_kind)%array(iparticle_local)
148 my_vel(:, iparticle) = particles%els(iparticle)%v(:)
153 CALL shake_control(gci=gci, local_molecules=local_molecules, &
154 molecule_set=molecules%els, molecule_kind_set=molecule_kinds%els, &
155 particle_set=particles%els, pos=my_pos, vel=my_vel, dt=mydt, &
156 shake_tol=shake_tol, log_unit=my_log_unit, lagrange_mult=my_lagrange_mult, &
157 dump_lm=my_dump_lm, cell=cell, group=force_env%para_env, &
158 local_particles=local_particles)
161 IF (
PRESENT(reset))
THEN
164 DO i = 1,
SIZE(molecules%els)
165 IF (
ASSOCIATED(molecules%els(i)%lci%lcolv))
THEN
166 DO j = 1,
SIZE(molecules%els(i)%lci%lcolv)
168 molecules%els(i)%lci%lcolv(j)%lambda = 0.0_dp
171 IF (
ASSOCIATED(molecules%els(i)%lci%lg3x3))
THEN
172 DO j = 1,
SIZE(molecules%els(i)%lci%lg3x3)
174 molecules%els(i)%lci%lg3x3(j)%lambda = 0.0_dp
177 IF (
ASSOCIATED(molecules%els(i)%lci%lg4x6))
THEN
178 DO j = 1,
SIZE(molecules%els(i)%lci%lg4x6)
180 molecules%els(i)%lci%lg4x6(j)%lambda = 0.0_dp
185 IF (
ASSOCIATED(gci))
THEN
186 IF (
ASSOCIATED(gci%lcolv))
THEN
187 DO j = 1,
SIZE(gci%lcolv)
189 gci%lcolv(j)%lambda = 0.0_dp
192 IF (
ASSOCIATED(gci%lg3x3))
THEN
193 DO j = 1,
SIZE(gci%lg3x3)
195 gci%lg3x3(j)%lambda = 0.0_dp
198 IF (
ASSOCIATED(gci%lg4x6))
THEN
199 DO j = 1,
SIZE(gci%lg4x6)
201 gci%lg4x6(j)%lambda = 0.0_dp
216 CALL timestop(handle)
238 TYPE(force_env_type),
POINTER :: force_env
239 REAL(kind=
dp),
INTENT(in),
OPTIONAL :: dt
240 REAL(kind=
dp),
INTENT(in) :: shake_tol
241 INTEGER,
INTENT(in),
OPTIONAL :: log_unit, lagrange_mult
242 LOGICAL,
INTENT(IN),
OPTIONAL :: dump_lm
243 REAL(kind=
dp),
DIMENSION(:, :),
INTENT(INOUT), &
244 OPTIONAL,
TARGET :: vel
245 LOGICAL,
INTENT(IN),
OPTIONAL :: reset
247 CHARACTER(len=*),
PARAMETER :: routinen =
'force_env_rattle'
249 INTEGER :: handle, i, iparticle, iparticle_kind, iparticle_local, j, my_lagrange_mult, &
250 my_log_unit, nparticle_kind, nparticle_local
251 LOGICAL :: has_vel, my_dump_lm
252 REAL(kind=
dp) :: mydt
253 REAL(kind=
dp),
DIMENSION(:, :),
POINTER :: my_vel
254 TYPE(atomic_kind_list_type),
POINTER :: atomic_kinds
255 TYPE(cell_type),
POINTER :: cell
256 TYPE(cp_subsys_type),
POINTER :: subsys
257 TYPE(distribution_1d_type),
POINTER :: local_molecules, local_particles
258 TYPE(global_constraint_type),
POINTER :: gci
259 TYPE(molecule_kind_list_type),
POINTER :: molecule_kinds
260 TYPE(molecule_list_type),
POINTER :: molecules
261 TYPE(particle_list_type),
POINTER :: particles
263 CALL timeset(routinen, handle)
264 cpassert(
ASSOCIATED(force_env))
265 cpassert(force_env%ref_count > 0)
267 IF (
PRESENT(log_unit)) my_log_unit = log_unit
268 my_lagrange_mult = -1
269 IF (
PRESENT(lagrange_mult)) my_lagrange_mult = lagrange_mult
271 IF (
PRESENT(dump_lm)) my_dump_lm = dump_lm
272 NULLIFY (subsys, cell, molecules, molecule_kinds, local_molecules, particles, &
274 IF (
PRESENT(vel)) my_vel => vel
276 IF (
PRESENT(dt)) mydt = dt
279 atomic_kinds=atomic_kinds, &
280 local_molecules=local_molecules, &
281 local_particles=local_particles, &
282 molecules=molecules, &
283 molecule_kinds=molecule_kinds, &
284 particles=particles, &
286 nparticle_kind = atomic_kinds%n_els
288 IF (.NOT.
ASSOCIATED(my_vel))
THEN
290 ALLOCATE (my_vel(3, particles%n_els))
292 DO iparticle_kind = 1, nparticle_kind
293 nparticle_local = local_particles%n_el(iparticle_kind)
294 DO iparticle_local = 1, nparticle_local
295 iparticle = local_particles%list(iparticle_kind)%array(iparticle_local)
296 my_vel(:, iparticle) = particles%els(iparticle)%v(:)
302 molecule_set=molecules%els, molecule_kind_set=molecule_kinds%els, &
303 particle_set=particles%els, vel=my_vel, dt=mydt, &
304 rattle_tol=shake_tol, log_unit=my_log_unit, lagrange_mult=my_lagrange_mult, &
305 dump_lm=my_dump_lm, cell=cell, group=force_env%para_env, &
306 local_particles=local_particles)
309 IF (
PRESENT(reset))
THEN
312 DO i = 1,
SIZE(molecules%els)
313 IF (
ASSOCIATED(molecules%els(i)%lci%lcolv))
THEN
314 DO j = 1,
SIZE(molecules%els(i)%lci%lcolv)
316 molecules%els(i)%lci%lcolv(j)%lambda = 0.0_dp
319 IF (
ASSOCIATED(molecules%els(i)%lci%lg3x3))
THEN
320 DO j = 1,
SIZE(molecules%els(i)%lci%lg3x3)
322 molecules%els(i)%lci%lg3x3(j)%lambda = 0.0_dp
325 IF (
ASSOCIATED(molecules%els(i)%lci%lg4x6))
THEN
326 DO j = 1,
SIZE(molecules%els(i)%lci%lg4x6)
328 molecules%els(i)%lci%lg4x6(j)%lambda = 0.0_dp
333 IF (
ASSOCIATED(gci))
THEN
334 IF (
ASSOCIATED(gci%lcolv))
THEN
335 DO j = 1,
SIZE(gci%lcolv)
337 gci%lcolv(j)%lambda = 0.0_dp
340 IF (
ASSOCIATED(gci%lg3x3))
THEN
341 DO j = 1,
SIZE(gci%lg3x3)
343 gci%lg3x3(j)%lambda = 0.0_dp
346 IF (
ASSOCIATED(gci%lg4x6))
THEN
347 DO j = 1,
SIZE(gci%lg4x6)
349 gci%lg4x6(j)%lambda = 0.0_dp
360 CALL timestop(handle)
369 TYPE(force_env_type),
POINTER :: force_env
371 CHARACTER(len=*),
PARAMETER :: routinen =
'rescale_forces'
373 INTEGER :: handle, iparticle
375 REAL(kind=
dp) :: force(3), max_value, mod_force
376 TYPE(cp_subsys_type),
POINTER :: subsys
377 TYPE(particle_list_type),
POINTER :: particles
378 TYPE(section_vals_type),
POINTER :: rescale_force_section
380 CALL timeset(routinen, handle)
381 cpassert(
ASSOCIATED(force_env))
382 cpassert(force_env%ref_count > 0)
389 DO iparticle = 1,
SIZE(particles%els)
390 force = particles%els(iparticle)%f(:)
391 mod_force = sqrt(dot_product(force, force))
392 IF ((mod_force > max_value) .AND. (mod_force /= 0.0_dp))
THEN
393 force = force/mod_force*max_value
394 particles%els(iparticle)%f(:) = force
398 CALL timestop(handle)
414 grand_total_force, zero_force_core_shell_atom)
416 TYPE(particle_list_type),
POINTER :: particles
417 INTEGER,
INTENT(IN) :: iw
418 CHARACTER(LEN=*),
INTENT(IN) :: label
419 INTEGER,
INTENT(IN) :: ndigits
420 REAL(kind=
dp),
DIMENSION(3),
INTENT(OUT) :: total_force
421 REAL(kind=
dp),
DIMENSION(3),
INTENT(INOUT), &
422 OPTIONAL :: grand_total_force
423 LOGICAL,
INTENT(IN),
OPTIONAL :: zero_force_core_shell_atom
425 CHARACTER(LEN=23) :: fmtstr3
426 CHARACTER(LEN=36) :: fmtstr2
427 CHARACTER(LEN=46) :: fmtstr1
428 INTEGER :: i, ikind, iparticle, n
429 LOGICAL :: zero_force
430 REAL(kind=
dp),
DIMENSION(3) :: f
433 cpassert(
ASSOCIATED(particles))
434 n = min(max(1, ndigits), 20)
435 fmtstr1 =
"(/,T2,A,/,/,T2,A,T11,A,T18,A,T35,A1,2( X,A1))"
436 WRITE (unit=fmtstr1(39:40), fmt=
"(I2)") n + 6
437 fmtstr2 =
"(T2,I6,1X,I6,T21,A,T28,3(1X,F . ))"
438 WRITE (unit=fmtstr2(33:34), fmt=
"(I2)") n
439 WRITE (unit=fmtstr2(30:31), fmt=
"(I2)") n + 6
440 fmtstr3 =
"(T2,A,T28,4(1X,F . ))"
441 WRITE (unit=fmtstr3(20:21), fmt=
"(I2)") n
442 WRITE (unit=fmtstr3(17:18), fmt=
"(I2)") n + 6
443 IF (
PRESENT(zero_force_core_shell_atom))
THEN
444 zero_force = zero_force_core_shell_atom
448 WRITE (unit=iw, fmt=fmtstr1) &
449 label//
" FORCES in [a.u.]",
"# Atom",
"Kind",
"Element",
"X",
"Y",
"Z"
450 total_force(1:3) = 0.0_dp
451 DO iparticle = 1, particles%n_els
452 ikind = particles%els(iparticle)%atomic_kind%kind_number
453 IF (particles%els(iparticle)%atom_index /= 0)
THEN
454 i = particles%els(iparticle)%atom_index
458 IF (zero_force .AND. (particles%els(iparticle)%shell_index /= 0))
THEN
461 f(1:3) = particles%els(iparticle)%f(1:3)
463 WRITE (unit=iw, fmt=fmtstr2) &
464 i, ikind, particles%els(iparticle)%atomic_kind%element_symbol, f(1:3)
465 total_force(1:3) = total_force(1:3) + f(1:3)
467 WRITE (unit=iw, fmt=fmtstr3) &
468 "SUM OF "//label//
" FORCES", total_force(1:3), sqrt(sum(total_force(:)**2))
471 IF (
PRESENT(grand_total_force))
THEN
472 grand_total_force(1:3) = grand_total_force(1:3) + total_force(1:3)
473 WRITE (unit=iw, fmt=
"(A)")
""
474 WRITE (unit=iw, fmt=fmtstr3) &
475 "GRAND TOTAL FORCE", grand_total_force(1:3), sqrt(sum(grand_total_force(:)**2))
492 INTEGER,
INTENT(IN) :: iounit
493 TYPE(particle_list_type),
POINTER :: particles
494 REAL(kind=
dp),
DIMENSION(:),
INTENT(IN) :: atener
495 CHARACTER(LEN=*),
INTENT(IN) :: label
500 WRITE (unit=iounit, fmt=
"(/,T2,A)") trim(label)
501 WRITE (unit=iounit, fmt=
"(T4,A,T30,A,T42,A,T54,A,T69,A)") &
502 "Atom Element",
"X",
"Y",
"Z",
"Energy[a.u.]"
503 natom = particles%n_els
505 WRITE (iounit,
"(I6,T12,A2,T24,3F12.6,F21.10)") i, &
506 trim(adjustl(particles%els(i)%atomic_kind%element_symbol)), &
507 particles%els(i)%r(1:3)*
angstrom, atener(i)
509 WRITE (iounit,
"(A)")
""
represent a simple array based list of the given type
Handles all functions related to the CELL.
Contains routines useful for the application of constraints during MD.
subroutine, public getold(gci, local_molecules, molecule_set, molecule_kind_set, particle_set, cell)
saves all of the old variables
subroutine, public rattle_control(gci, local_molecules, molecule_set, molecule_kind_set, particle_set, vel, dt, rattle_tol, log_unit, lagrange_mult, dump_lm, cell, group, local_particles)
...
subroutine, public shake_control(gci, local_molecules, molecule_set, molecule_kind_set, particle_set, pos, vel, dt, shake_tol, log_unit, lagrange_mult, dump_lm, cell, group, local_particles)
...
types that represent a subsys, i.e. a part of the system
subroutine, public cp_subsys_get(subsys, ref_count, atomic_kinds, atomic_kind_set, particles, particle_set, local_particles, molecules, molecule_set, molecule_kinds, molecule_kind_set, local_molecules, para_env, colvar_p, shell_particles, core_particles, gci, multipoles, natom, nparticle, ncore, nshell, nkind, atprop, virial, results, cell)
returns information about various attributes of the given subsys
stores a lists of integer that are local to a processor. The idea is that these integers represent ob...
Interface for the force calculations.
recursive subroutine, public force_env_get(force_env, in_use, fist_env, qs_env, meta_env, fp_env, subsys, para_env, potential_energy, additional_potential, kinetic_energy, harmonic_shell, kinetic_shell, cell, sub_force_env, qmmm_env, qmmmx_env, eip_env, pwdft_env, globenv, input, force_env_section, method_name_id, root_section, mixed_env, nnp_env, embed_env)
returns various attributes about the force environment
subroutine, public write_forces(particles, iw, label, ndigits, total_force, grand_total_force, zero_force_core_shell_atom)
Write forces.
subroutine, public force_env_shake(force_env, dt, shake_tol, log_unit, lagrange_mult, dump_lm, pos, vel, compold, reset)
perform shake (enforcing of constraints)
subroutine, public write_atener(iounit, particles, atener, label)
Write the atomic coordinates, types, and energies.
subroutine, public rescale_forces(force_env)
Rescale forces if requested.
subroutine, public force_env_rattle(force_env, dt, shake_tol, log_unit, lagrange_mult, dump_lm, vel, reset)
perform rattle (enforcing of constraints on velocities) This routine can be easily adapted to perform...
Defines the basic variable types.
integer, parameter, public dp
represent a simple array based list of the given type
represent a simple array based list of the given type
Define the data structure for the molecule information.
represent a simple array based list of the given type
Define the data structure for the particle information.
subroutine, public update_particle_set(particle_set, int_group, pos, vel, for, add)
...
Definition of physical constants:
real(kind=dp), parameter, public angstrom