-
Notifications
You must be signed in to change notification settings - Fork 1
/
qs_force.F
683 lines (619 loc) · 34.6 KB
/
qs_force.F
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
509
510
511
512
513
514
515
516
517
518
519
520
521
522
523
524
525
526
527
528
529
530
531
532
533
534
535
536
537
538
539
540
541
542
543
544
545
546
547
548
549
550
551
552
553
554
555
556
557
558
559
560
561
562
563
564
565
566
567
568
569
570
571
572
573
574
575
576
577
578
579
580
581
582
583
584
585
586
587
588
589
590
591
592
593
594
595
596
597
598
599
600
601
602
603
604
605
606
607
608
609
610
611
612
613
614
615
616
617
618
619
620
621
622
623
624
625
626
627
628
629
630
631
632
633
634
635
636
637
638
639
640
641
642
643
644
645
646
647
648
649
650
651
652
653
654
655
656
657
658
659
660
661
662
663
664
665
666
667
668
669
670
671
672
673
674
675
676
677
678
679
680
681
682
683
!--------------------------------------------------------------------------------------------------!
! CP2K: A general program to perform molecular dynamics simulations !
! Copyright 2000-2024 CP2K developers group <https://cp2k.org> !
! !
! SPDX-License-Identifier: GPL-2.0-or-later !
!--------------------------------------------------------------------------------------------------!
! **************************************************************************************************
!> \brief Quickstep force driver routine
!> \author MK (12.06.2002)
! **************************************************************************************************
MODULE qs_force
USE atomic_kind_types, ONLY: atomic_kind_type,&
get_atomic_kind_set
USE cp_control_types, ONLY: dft_control_type
USE cp_dbcsr_api, ONLY: dbcsr_copy,&
dbcsr_p_type,&
dbcsr_set
USE cp_dbcsr_operations, ONLY: dbcsr_allocate_matrix_set,&
dbcsr_deallocate_matrix_set
USE cp_dbcsr_output, ONLY: cp_dbcsr_write_sparse_matrix
USE cp_log_handling, ONLY: cp_get_default_logger,&
cp_logger_get_default_io_unit,&
cp_logger_type
USE cp_output_handling, ONLY: cp_p_file,&
cp_print_key_finished_output,&
cp_print_key_should_output,&
cp_print_key_unit_nr
USE dft_plus_u, ONLY: plus_u
USE ec_env_types, ONLY: energy_correction_type
USE efield_utils, ONLY: calculate_ecore_efield,&
efield_potential_lengh_gauge
USE energy_corrections, ONLY: energy_correction
USE excited_states, ONLY: excited_state_energy
USE hfx_exx, ONLY: calculate_exx
USE input_constants, ONLY: ri_mp2_laplace,&
ri_mp2_method_gpw,&
ri_rpa_method_gpw
USE input_section_types, ONLY: section_vals_get,&
section_vals_get_subs_vals,&
section_vals_type,&
section_vals_val_get
USE kinds, ONLY: dp
USE lri_environment_types, ONLY: lri_environment_type
USE message_passing, ONLY: mp_para_env_type
USE mp2_cphf, ONLY: update_mp2_forces
USE mulliken, ONLY: mulliken_restraint
USE particle_types, ONLY: particle_type
USE qs_core_energies, ONLY: calculate_ecore_overlap,&
calculate_ecore_self
USE qs_core_hamiltonian, ONLY: build_core_hamiltonian_matrix
USE qs_dftb_dispersion, ONLY: calculate_dftb_dispersion
USE qs_dftb_matrices, ONLY: build_dftb_matrices
USE qs_energy, ONLY: qs_energies
USE qs_energy_types, ONLY: qs_energy_type
USE qs_environment_methods, ONLY: qs_env_rebuild_pw_env
USE qs_environment_types, ONLY: get_qs_env,&
qs_environment_type
USE qs_external_potential, ONLY: external_c_potential,&
external_e_potential
USE qs_force_types, ONLY: allocate_qs_force,&
qs_force_type,&
replicate_qs_force,&
zero_qs_force
USE qs_ks_methods, ONLY: qs_ks_update_qs_env
USE qs_ks_types, ONLY: qs_ks_env_type,&
set_ks_env
USE qs_rho_types, ONLY: qs_rho_get,&
qs_rho_type
USE qs_scf_post_scf, ONLY: qs_scf_compute_properties
USE qs_subsys_types, ONLY: qs_subsys_set,&
qs_subsys_type
USE ri_environment_methods, ONLY: build_ri_matrices
USE rt_propagation_forces, ONLY: calc_c_mat_force,&
rt_admm_force
USE rt_propagation_velocity_gauge, ONLY: velocity_gauge_ks_matrix,&
velocity_gauge_nl_force
USE se_core_core, ONLY: se_core_core_interaction
USE se_core_matrix, ONLY: build_se_core_matrix
USE virial_types, ONLY: symmetrize_virial,&
virial_type
USE xtb_matrices, ONLY: build_xtb_matrices
#include "./base/base_uses.f90"
IMPLICIT NONE
PRIVATE
! *** Global parameters ***
CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'qs_force'
! *** Public subroutines ***
PUBLIC :: qs_calc_energy_force
CONTAINS
! **************************************************************************************************
!> \brief ...
!> \param qs_env ...
!> \param calc_force ...
!> \param consistent_energies ...
!> \param linres ...
! **************************************************************************************************
SUBROUTINE qs_calc_energy_force(qs_env, calc_force, consistent_energies, linres)
TYPE(qs_environment_type), POINTER :: qs_env
LOGICAL :: calc_force, consistent_energies, linres
qs_env%linres_run = linres
IF (calc_force) THEN
CALL qs_forces(qs_env)
ELSE
CALL qs_energies(qs_env, calc_forces=.FALSE., &
consistent_energies=consistent_energies)
END IF
END SUBROUTINE qs_calc_energy_force
! **************************************************************************************************
!> \brief Calculate the Quickstep forces.
!> \param qs_env ...
!> \date 29.10.2002
!> \author MK
!> \version 1.0
! **************************************************************************************************
SUBROUTINE qs_forces(qs_env)
TYPE(qs_environment_type), POINTER :: qs_env
CHARACTER(len=*), PARAMETER :: routineN = 'qs_forces'
INTEGER :: after, handle, i, iatom, ic, ikind, &
ispin, iw, natom, nkind, nspin, &
output_unit
INTEGER, ALLOCATABLE, DIMENSION(:) :: atom_of_kind, kind_of, natom_of_kind
LOGICAL :: do_admm, do_exx, do_gw, do_im_time, &
has_unit_metric, omit_headers, &
perform_ec, reuse_hfx
REAL(dp) :: dummy_real, dummy_real2(2)
TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
TYPE(cp_logger_type), POINTER :: logger
TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_s, matrix_w, rho_ao
TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: matrix_w_kp
TYPE(dft_control_type), POINTER :: dft_control
TYPE(energy_correction_type), POINTER :: ec_env
TYPE(lri_environment_type), POINTER :: lri_env
TYPE(mp_para_env_type), POINTER :: para_env
TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
TYPE(qs_energy_type), POINTER :: energy
TYPE(qs_force_type), DIMENSION(:), POINTER :: force
TYPE(qs_ks_env_type), POINTER :: ks_env
TYPE(qs_rho_type), POINTER :: rho
TYPE(qs_subsys_type), POINTER :: subsys
TYPE(section_vals_type), POINTER :: hfx_sections, print_section
TYPE(virial_type), POINTER :: virial
CALL timeset(routineN, handle)
NULLIFY (logger)
logger => cp_get_default_logger()
! rebuild plane wave environment
CALL qs_env_rebuild_pw_env(qs_env)
! zero out the forces in particle set
CALL get_qs_env(qs_env, particle_set=particle_set)
natom = SIZE(particle_set)
DO iatom = 1, natom
particle_set(iatom)%f = 0.0_dp
END DO
! get atom mapping
NULLIFY (atomic_kind_set)
CALL get_qs_env(qs_env, atomic_kind_set=atomic_kind_set)
CALL get_atomic_kind_set(atomic_kind_set=atomic_kind_set, &
atom_of_kind=atom_of_kind, &
kind_of=kind_of)
NULLIFY (force, subsys, dft_control)
CALL get_qs_env(qs_env, &
force=force, &
subsys=subsys, &
dft_control=dft_control)
IF (.NOT. ASSOCIATED(force)) THEN
! *** Allocate the force data structure ***
nkind = SIZE(atomic_kind_set)
CALL get_atomic_kind_set(atomic_kind_set=atomic_kind_set, natom_of_kind=natom_of_kind)
CALL allocate_qs_force(force, natom_of_kind)
DEALLOCATE (natom_of_kind)
CALL qs_subsys_set(subsys, force=force)
END IF
CALL zero_qs_force(force)
! Check if CDFT potential is needed and save it until forces have been calculated
IF (dft_control%qs_control%cdft) THEN
dft_control%qs_control%cdft_control%save_pot = .TRUE.
END IF
! recalculate energy with forces
CALL qs_energies(qs_env, calc_forces=.TRUE.)
NULLIFY (para_env)
CALL get_qs_env(qs_env, &
para_env=para_env)
! Now we handle some special cases
! Maybe some of these would be better dealt with in qs_energies?
IF (qs_env%run_rtp) THEN
NULLIFY (matrix_w, matrix_s, ks_env)
CALL get_qs_env(qs_env, &
ks_env=ks_env, &
matrix_w=matrix_w, &
matrix_s=matrix_s)
CALL dbcsr_allocate_matrix_set(matrix_w, dft_control%nspins)
DO ispin = 1, dft_control%nspins
ALLOCATE (matrix_w(ispin)%matrix)
CALL dbcsr_copy(matrix_w(ispin)%matrix, matrix_s(1)%matrix, &
name="W MATRIX")
CALL dbcsr_set(matrix_w(ispin)%matrix, 0.0_dp)
END DO
CALL set_ks_env(ks_env, matrix_w=matrix_w)
CALL calc_c_mat_force(qs_env)
IF (dft_control%do_admm) CALL rt_admm_force(qs_env)
IF (dft_control%rtp_control%velocity_gauge .AND. dft_control%rtp_control%nl_gauge_transform) &
CALL velocity_gauge_nl_force(qs_env, particle_set)
END IF
! from an eventual Mulliken restraint
IF (dft_control%qs_control%mulliken_restraint) THEN
NULLIFY (matrix_w, matrix_s, rho)
CALL get_qs_env(qs_env, &
matrix_w=matrix_w, &
matrix_s=matrix_s, &
rho=rho)
NULLIFY (rho_ao)
CALL qs_rho_get(rho, rho_ao=rho_ao)
CALL mulliken_restraint(dft_control%qs_control%mulliken_restraint_control, &
para_env, matrix_s(1)%matrix, rho_ao, w_matrix=matrix_w)
END IF
! Add non-Pulay contribution of DFT+U to W matrix, since it has also to be
! digested with overlap matrix derivatives
IF (dft_control%dft_plus_u) THEN
NULLIFY (matrix_w)
CALL get_qs_env(qs_env, matrix_w=matrix_w)
CALL plus_u(qs_env=qs_env, matrix_w=matrix_w)
END IF
! Write W Matrix to output (if requested)
CALL get_qs_env(qs_env, has_unit_metric=has_unit_metric)
IF (.NOT. has_unit_metric) THEN
NULLIFY (matrix_w_kp)
CALL get_qs_env(qs_env, matrix_w_kp=matrix_w_kp)
nspin = SIZE(matrix_w_kp, 1)
DO ispin = 1, nspin
IF (BTEST(cp_print_key_should_output(logger%iter_info, &
qs_env%input, "DFT%PRINT%AO_MATRICES/W_MATRIX"), cp_p_file)) THEN
iw = cp_print_key_unit_nr(logger, qs_env%input, "DFT%PRINT%AO_MATRICES/W_MATRIX", &
extension=".Log")
CALL section_vals_val_get(qs_env%input, "DFT%PRINT%AO_MATRICES%NDIGITS", i_val=after)
CALL section_vals_val_get(qs_env%input, "DFT%PRINT%AO_MATRICES%OMIT_HEADERS", l_val=omit_headers)
after = MIN(MAX(after, 1), 16)
DO ic = 1, SIZE(matrix_w_kp, 2)
CALL cp_dbcsr_write_sparse_matrix(matrix_w_kp(ispin, ic)%matrix, 4, after, qs_env, &
para_env, output_unit=iw, omit_headers=omit_headers)
END DO
CALL cp_print_key_finished_output(iw, logger, qs_env%input, &
"DFT%PRINT%AO_MATRICES/W_MATRIX")
END IF
END DO
END IF
! Check if energy correction should be skipped
perform_ec = .FALSE.
IF (qs_env%energy_correction) THEN
CALL get_qs_env(qs_env, ec_env=ec_env)
IF (.NOT. ec_env%do_skip) perform_ec = .TRUE.
END IF
! Compute core forces (also overwrites matrix_w)
IF (dft_control%qs_control%semi_empirical) THEN
CALL build_se_core_matrix(qs_env=qs_env, para_env=para_env, &
calculate_forces=.TRUE.)
CALL se_core_core_interaction(qs_env, para_env, calculate_forces=.TRUE.)
ELSEIF (dft_control%qs_control%dftb) THEN
CALL build_dftb_matrices(qs_env=qs_env, para_env=para_env, &
calculate_forces=.TRUE.)
CALL calculate_dftb_dispersion(qs_env=qs_env, para_env=para_env, &
calculate_forces=.TRUE.)
ELSEIF (dft_control%qs_control%xtb) THEN
CALL build_xtb_matrices(qs_env=qs_env, calculate_forces=.TRUE.)
ELSEIF (perform_ec) THEN
! Calculates core and grid based forces
CALL energy_correction(qs_env, ec_init=.FALSE., calculate_forces=.TRUE.)
ELSE
! Dispersion energy and forces are calculated in qs_energy?
CALL build_core_hamiltonian_matrix(qs_env=qs_env, calculate_forces=.TRUE.)
! The above line reset the core H, which should be re-updated in case a TD field is applied:
IF (qs_env%run_rtp) THEN
IF (dft_control%apply_efield_field) &
CALL efield_potential_lengh_gauge(qs_env)
IF (dft_control%rtp_control%velocity_gauge) &
CALL velocity_gauge_ks_matrix(qs_env, subtract_nl_term=.FALSE.)
END IF
CALL calculate_ecore_self(qs_env)
CALL calculate_ecore_overlap(qs_env, para_env, calculate_forces=.TRUE.)
CALL calculate_ecore_efield(qs_env, calculate_forces=.TRUE.)
!swap external_e_potential before external_c_potential, to ensure
!that external potential on grid is loaded before calculating energy of cores
CALL external_e_potential(qs_env)
IF (.NOT. dft_control%qs_control%gapw) THEN
CALL external_c_potential(qs_env, calculate_forces=.TRUE.)
END IF
! RIGPW matrices
IF (dft_control%qs_control%rigpw) THEN
CALL get_qs_env(qs_env=qs_env, lri_env=lri_env)
CALL build_ri_matrices(lri_env, qs_env, calculate_forces=.TRUE.)
END IF
END IF
! MP2 Code
IF (ASSOCIATED(qs_env%mp2_env)) THEN
NULLIFY (energy)
CALL get_qs_env(qs_env, energy=energy)
CALL qs_scf_compute_properties(qs_env, wf_type='MP2 ', do_mp2=.TRUE.)
CALL qs_ks_update_qs_env(qs_env, just_energy=.TRUE.)
energy%total = energy%total + energy%mp2
IF ((qs_env%mp2_env%method == ri_mp2_method_gpw .OR. qs_env%mp2_env%method == ri_mp2_laplace .OR. &
qs_env%mp2_env%method == ri_rpa_method_gpw) &
.AND. .NOT. qs_env%mp2_env%do_im_time) THEN
CALL update_mp2_forces(qs_env)
END IF
!RPA EXX energy and forces
IF (qs_env%mp2_env%method == ri_rpa_method_gpw) THEN
do_exx = .FALSE.
hfx_sections => section_vals_get_subs_vals(qs_env%input, "DFT%XC%WF_CORRELATION%RI_RPA%HF")
CALL section_vals_get(hfx_sections, explicit=do_exx)
IF (do_exx) THEN
do_gw = qs_env%mp2_env%ri_rpa%do_ri_g0w0
do_admm = qs_env%mp2_env%ri_rpa%do_admm
reuse_hfx = qs_env%mp2_env%ri_rpa%reuse_hfx
do_im_time = qs_env%mp2_env%do_im_time
output_unit = cp_logger_get_default_io_unit()
dummy_real = 0.0_dp
CALL calculate_exx(qs_env=qs_env, &
unit_nr=output_unit, &
hfx_sections=hfx_sections, &
x_data=qs_env%mp2_env%ri_rpa%x_data, &
do_gw=do_gw, &
do_admm=do_admm, &
calc_forces=.TRUE., &
reuse_hfx=reuse_hfx, &
do_im_time=do_im_time, &
E_ex_from_GW=dummy_real, &
E_admm_from_GW=dummy_real2, &
t3=dummy_real)
END IF
END IF
ELSEIF (perform_ec) THEN
! energy correction forces postponed
ELSEIF (qs_env%harris_method) THEN
! Harris method forces already done in harris_energy_correction
ELSE
! Compute grid-based forces
CALL qs_ks_update_qs_env(qs_env, calculate_forces=.TRUE.)
END IF
! Excited state forces
CALL excited_state_energy(qs_env, calculate_forces=.TRUE.)
! replicate forces (get current pointer)
NULLIFY (force)
CALL get_qs_env(qs_env=qs_env, force=force)
CALL replicate_qs_force(force, para_env)
DO iatom = 1, natom
ikind = kind_of(iatom)
i = atom_of_kind(iatom)
! XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX
! the force is - dE/dR, what is called force is actually the gradient
! Things should have the right name
! The minus sign below is a hack
! XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX
force(ikind)%other(1:3, i) = -particle_set(iatom)%f(1:3) + force(ikind)%ch_pulay(1:3, i)
force(ikind)%total(1:3, i) = force(ikind)%total(1:3, i) + force(ikind)%other(1:3, i)
particle_set(iatom)%f = -force(ikind)%total(1:3, i)
END DO
NULLIFY (virial, energy)
CALL get_qs_env(qs_env=qs_env, virial=virial, energy=energy)
IF (virial%pv_availability) THEN
CALL para_env%sum(virial%pv_overlap)
CALL para_env%sum(virial%pv_ekinetic)
CALL para_env%sum(virial%pv_ppl)
CALL para_env%sum(virial%pv_ppnl)
CALL para_env%sum(virial%pv_ecore_overlap)
CALL para_env%sum(virial%pv_ehartree)
CALL para_env%sum(virial%pv_exc)
CALL para_env%sum(virial%pv_exx)
CALL para_env%sum(virial%pv_vdw)
CALL para_env%sum(virial%pv_mp2)
CALL para_env%sum(virial%pv_nlcc)
CALL para_env%sum(virial%pv_gapw)
CALL para_env%sum(virial%pv_lrigpw)
CALL para_env%sum(virial%pv_virial)
CALL symmetrize_virial(virial)
! Add the volume terms of the virial
IF ((.NOT. virial%pv_numer) .AND. &
(.NOT. (dft_control%qs_control%dftb .OR. &
dft_control%qs_control%xtb .OR. &
dft_control%qs_control%semi_empirical))) THEN
! Harris energy correction requires volume terms from
! 1) Harris functional contribution, and
! 2) Linear Response solver
IF (perform_ec) THEN
CALL get_qs_env(qs_env, ec_env=ec_env)
energy%hartree = ec_env%ehartree
energy%exc = ec_env%exc
IF (dft_control%do_admm) THEN
energy%exc_aux_fit = ec_env%exc_aux_fit
END IF
END IF
DO i = 1, 3
virial%pv_ehartree(i, i) = virial%pv_ehartree(i, i) &
- 2.0_dp*(energy%hartree + energy%sccs_pol)
virial%pv_virial(i, i) = virial%pv_virial(i, i) - energy%exc &
- 2.0_dp*(energy%hartree + energy%sccs_pol)
virial%pv_exc(i, i) = virial%pv_exc(i, i) - energy%exc
IF (dft_control%do_admm) THEN
virial%pv_exc(i, i) = virial%pv_exc(i, i) - energy%exc_aux_fit
virial%pv_virial(i, i) = virial%pv_virial(i, i) - energy%exc_aux_fit
END IF
! The factor 2 is a hack. It compensates the plus sign in h_stress/pw_poisson_solve.
! The sign in pw_poisson_solve is correct for FIST, but not for QS.
! There should be a more elegant solution to that ...
END DO
END IF
END IF
output_unit = cp_print_key_unit_nr(logger, qs_env%input, "DFT%PRINT%DERIVATIVES", &
extension=".Log")
print_section => section_vals_get_subs_vals(qs_env%input, "DFT%PRINT%DERIVATIVES")
IF (dft_control%qs_control%semi_empirical) THEN
CALL write_forces(force, atomic_kind_set, 2, output_unit=output_unit, &
print_section=print_section)
ELSE IF (dft_control%qs_control%dftb) THEN
CALL write_forces(force, atomic_kind_set, 4, output_unit=output_unit, &
print_section=print_section)
ELSE IF (dft_control%qs_control%xtb) THEN
CALL write_forces(force, atomic_kind_set, 4, output_unit=output_unit, &
print_section=print_section)
ELSE IF (dft_control%qs_control%gapw .OR. dft_control%qs_control%gapw_xc) THEN
CALL write_forces(force, atomic_kind_set, 1, output_unit=output_unit, &
print_section=print_section)
ELSE
CALL write_forces(force, atomic_kind_set, 0, output_unit=output_unit, &
print_section=print_section)
END IF
CALL cp_print_key_finished_output(output_unit, logger, qs_env%input, &
"DFT%PRINT%DERIVATIVES")
! deallocate W Matrix:
NULLIFY (ks_env, matrix_w_kp)
CALL get_qs_env(qs_env=qs_env, &
matrix_w_kp=matrix_w_kp, &
ks_env=ks_env)
CALL dbcsr_deallocate_matrix_set(matrix_w_kp)
NULLIFY (matrix_w_kp)
CALL set_ks_env(ks_env, matrix_w_kp=matrix_w_kp)
DEALLOCATE (atom_of_kind, kind_of)
CALL timestop(handle)
END SUBROUTINE qs_forces
! **************************************************************************************************
!> \brief Write a Quickstep force data structure to output unit
!> \param qs_force ...
!> \param atomic_kind_set ...
!> \param ftype ...
!> \param output_unit ...
!> \param print_section ...
!> \date 05.06.2002
!> \author MK
!> \version 1.0
! **************************************************************************************************
SUBROUTINE write_forces(qs_force, atomic_kind_set, ftype, output_unit, &
print_section)
TYPE(qs_force_type), DIMENSION(:), POINTER :: qs_force
TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
INTEGER, INTENT(IN) :: ftype, output_unit
TYPE(section_vals_type), POINTER :: print_section
CHARACTER(LEN=13) :: fmtstr5
CHARACTER(LEN=15) :: fmtstr4
CHARACTER(LEN=20) :: fmtstr3
CHARACTER(LEN=35) :: fmtstr2
CHARACTER(LEN=48) :: fmtstr1
INTEGER :: i, iatom, ikind, my_ftype, natom, ndigits
INTEGER, ALLOCATABLE, DIMENSION(:) :: atom_of_kind, kind_of
REAL(KIND=dp), DIMENSION(3) :: grand_total
IF (output_unit > 0) THEN
IF (.NOT. ASSOCIATED(qs_force)) THEN
CALL cp_abort(__LOCATION__, &
"The qs_force pointer is not associated "// &
"and cannot be printed")
END IF
CALL get_atomic_kind_set(atomic_kind_set=atomic_kind_set, atom_of_kind=atom_of_kind, &
kind_of=kind_of, natom=natom)
! Variable precision output of the forces
CALL section_vals_val_get(print_section, "NDIGITS", &
i_val=ndigits)
fmtstr1 = "(/,/,T2,A,/,/,T3,A,T11,A,T23,A,T40,A1,2( X,A1))"
WRITE (UNIT=fmtstr1(41:42), FMT="(I2)") ndigits + 5
fmtstr2 = "(/,(T2,I5,4X,I4,T18,A,T34,3F . ))"
WRITE (UNIT=fmtstr2(32:33), FMT="(I2)") ndigits
WRITE (UNIT=fmtstr2(29:30), FMT="(I2)") ndigits + 6
fmtstr3 = "(/,T3,A,T34,3F . )"
WRITE (UNIT=fmtstr3(18:19), FMT="(I2)") ndigits
WRITE (UNIT=fmtstr3(15:16), FMT="(I2)") ndigits + 6
fmtstr4 = "((T34,3F . ))"
WRITE (UNIT=fmtstr4(12:13), FMT="(I2)") ndigits
WRITE (UNIT=fmtstr4(9:10), FMT="(I2)") ndigits + 6
fmtstr5 = "(/T2,A//T3,A)"
WRITE (UNIT=output_unit, FMT=fmtstr1) &
"FORCES [a.u.]", "Atom", "Kind", "Component", "X", "Y", "Z"
grand_total(:) = 0.0_dp
my_ftype = ftype
SELECT CASE (my_ftype)
CASE DEFAULT
DO iatom = 1, natom
ikind = kind_of(iatom)
i = atom_of_kind(iatom)
WRITE (UNIT=output_unit, FMT=fmtstr2) &
iatom, ikind, " total", qs_force(ikind)%total(1:3, i)
grand_total(1:3) = grand_total(1:3) + qs_force(ikind)%total(1:3, i)
END DO
CASE (0)
DO iatom = 1, natom
ikind = kind_of(iatom)
i = atom_of_kind(iatom)
WRITE (UNIT=output_unit, FMT=fmtstr2) &
iatom, ikind, " overlap", qs_force(ikind)%overlap(1:3, i), &
iatom, ikind, " overlap_admm", qs_force(ikind)%overlap_admm(1:3, i), &
iatom, ikind, " kinetic", qs_force(ikind)%kinetic(1:3, i), &
iatom, ikind, " gth_ppl", qs_force(ikind)%gth_ppl(1:3, i), &
iatom, ikind, " gth_nlcc", qs_force(ikind)%gth_nlcc(1:3, i), &
iatom, ikind, " gth_ppnl", qs_force(ikind)%gth_ppnl(1:3, i), &
iatom, ikind, " core_overlap", qs_force(ikind)%core_overlap(1:3, i), &
iatom, ikind, " rho_core", qs_force(ikind)%rho_core(1:3, i), &
iatom, ikind, " rho_elec", qs_force(ikind)%rho_elec(1:3, i), &
iatom, ikind, " rho_lri_elec", qs_force(ikind)%rho_lri_elec(1:3, i), &
iatom, ikind, " ch_pulay", qs_force(ikind)%ch_pulay(1:3, i), &
iatom, ikind, " dispersion", qs_force(ikind)%dispersion(1:3, i), &
iatom, ikind, " gCP", qs_force(ikind)%gcp(1:3, i), &
iatom, ikind, " other", qs_force(ikind)%other(1:3, i), &
iatom, ikind, " fock_4c", qs_force(ikind)%fock_4c(1:3, i), &
iatom, ikind, " ehrenfest", qs_force(ikind)%ehrenfest(1:3, i), &
iatom, ikind, " efield", qs_force(ikind)%efield(1:3, i), &
iatom, ikind, " eev", qs_force(ikind)%eev(1:3, i), &
iatom, ikind, " mp2_non_sep", qs_force(ikind)%mp2_non_sep(1:3, i), &
iatom, ikind, " total", qs_force(ikind)%total(1:3, i)
grand_total(1:3) = grand_total(1:3) + qs_force(ikind)%total(1:3, i)
END DO
CASE (1)
DO iatom = 1, natom
ikind = kind_of(iatom)
i = atom_of_kind(iatom)
WRITE (UNIT=output_unit, FMT=fmtstr2) &
iatom, ikind, " overlap", qs_force(ikind)%overlap(1:3, i), &
iatom, ikind, " overlap_admm", qs_force(ikind)%overlap_admm(1:3, i), &
iatom, ikind, " kinetic", qs_force(ikind)%kinetic(1:3, i), &
iatom, ikind, " gth_ppl", qs_force(ikind)%gth_ppl(1:3, i), &
iatom, ikind, " gth_nlcc", qs_force(ikind)%gth_nlcc(1:3, i), &
iatom, ikind, " gth_ppnl", qs_force(ikind)%gth_ppnl(1:3, i), &
iatom, ikind, " all_potential", qs_force(ikind)%all_potential(1:3, i), &
iatom, ikind, " core_overlap", qs_force(ikind)%core_overlap(1:3, i), &
iatom, ikind, " rho_core", qs_force(ikind)%rho_core(1:3, i), &
iatom, ikind, " rho_elec", qs_force(ikind)%rho_elec(1:3, i), &
iatom, ikind, " rho_lri_elec", qs_force(ikind)%rho_lri_elec(1:3, i), &
iatom, ikind, " vhxc_atom", qs_force(ikind)%vhxc_atom(1:3, i), &
iatom, ikind, " g0s_Vh_elec", qs_force(ikind)%g0s_Vh_elec(1:3, i), &
iatom, ikind, " ch_pulay", qs_force(ikind)%ch_pulay(1:3, i), &
iatom, ikind, " dispersion", qs_force(ikind)%dispersion(1:3, i), &
iatom, ikind, " gCP", qs_force(ikind)%gcp(1:3, i), &
iatom, ikind, " fock_4c", qs_force(ikind)%fock_4c(1:3, i), &
iatom, ikind, " ehrenfest", qs_force(ikind)%ehrenfest(1:3, i), &
iatom, ikind, " efield", qs_force(ikind)%efield(1:3, i), &
iatom, ikind, " eev", qs_force(ikind)%eev(1:3, i), &
iatom, ikind, " mp2_non_sep", qs_force(ikind)%mp2_non_sep(1:3, i), &
iatom, ikind, " total", qs_force(ikind)%total(1:3, i)
grand_total(1:3) = grand_total(1:3) + qs_force(ikind)%total(1:3, i)
END DO
CASE (2)
DO iatom = 1, natom
ikind = kind_of(iatom)
i = atom_of_kind(iatom)
WRITE (UNIT=output_unit, FMT=fmtstr2) &
iatom, ikind, " all_potential", qs_force(ikind)%all_potential(1:3, i), &
iatom, ikind, " rho_elec", qs_force(ikind)%rho_elec(1:3, i), &
iatom, ikind, " total", qs_force(ikind)%total(1:3, i)
grand_total(1:3) = grand_total(1:3) + qs_force(ikind)%total(1:3, i)
END DO
CASE (3)
DO iatom = 1, natom
ikind = kind_of(iatom)
i = atom_of_kind(iatom)
WRITE (UNIT=output_unit, FMT=fmtstr2) &
iatom, ikind, " overlap", qs_force(ikind)%overlap(1:3, i), &
iatom, ikind, "overlap_admm", qs_force(ikind)%overlap_admm(1:3, i), &
iatom, ikind, " kinetic", qs_force(ikind)%kinetic(1:3, i), &
iatom, ikind, " gth_ppl", qs_force(ikind)%gth_ppl(1:3, i), &
iatom, ikind, " gth_nlcc", qs_force(ikind)%gth_nlcc(1:3, i), &
iatom, ikind, " gth_ppnl", qs_force(ikind)%gth_ppnl(1:3, i), &
iatom, ikind, " core_overlap", qs_force(ikind)%core_overlap(1:3, i), &
iatom, ikind, " rho_core", qs_force(ikind)%rho_core(1:3, i), &
iatom, ikind, " rho_elec", qs_force(ikind)%rho_elec(1:3, i), &
iatom, ikind, " ch_pulay", qs_force(ikind)%ch_pulay(1:3, i), &
iatom, ikind, " fock_4c", qs_force(ikind)%fock_4c(1:3, i), &
iatom, ikind, " mp2_non_sep", qs_force(ikind)%mp2_non_sep(1:3, i), &
iatom, ikind, " total", qs_force(ikind)%total(1:3, i)
grand_total(1:3) = grand_total(1:3) + qs_force(ikind)%total(1:3, i)
END DO
CASE (4)
DO iatom = 1, natom
ikind = kind_of(iatom)
i = atom_of_kind(iatom)
WRITE (UNIT=output_unit, FMT=fmtstr2) &
iatom, ikind, " all_potential", qs_force(ikind)%all_potential(1:3, i), &
iatom, ikind, " overlap", qs_force(ikind)%overlap(1:3, i), &
iatom, ikind, " rho_elec", qs_force(ikind)%rho_elec(1:3, i), &
iatom, ikind, " repulsive", qs_force(ikind)%repulsive(1:3, i), &
iatom, ikind, " dispersion", qs_force(ikind)%dispersion(1:3, i), &
iatom, ikind, " efield", qs_force(ikind)%efield(1:3, i), &
iatom, ikind, " ehrenfest", qs_force(ikind)%ehrenfest(1:3, i), &
iatom, ikind, " total", qs_force(ikind)%total(1:3, i)
grand_total(1:3) = grand_total(1:3) + qs_force(ikind)%total(1:3, i)
END DO
CASE (5)
DO iatom = 1, natom
ikind = kind_of(iatom)
i = atom_of_kind(iatom)
WRITE (UNIT=output_unit, FMT=fmtstr2) &
iatom, ikind, " overlap", qs_force(ikind)%overlap(1:3, i), &
iatom, ikind, " kinetic", qs_force(ikind)%kinetic(1:3, i), &
iatom, ikind, " rho_elec", qs_force(ikind)%rho_elec(1:3, i), &
iatom, ikind, " dispersion", qs_force(ikind)%dispersion(1:3, i), &
iatom, ikind, " all potential", qs_force(ikind)%all_potential(1:3, i), &
iatom, ikind, " other", qs_force(ikind)%other(1:3, i), &
iatom, ikind, " total", qs_force(ikind)%total(1:3, i)
grand_total(1:3) = grand_total(1:3) + qs_force(ikind)%total(1:3, i)
END DO
END SELECT
WRITE (UNIT=output_unit, FMT=fmtstr3) "Sum of total", grand_total(1:3)
DEALLOCATE (atom_of_kind)
DEALLOCATE (kind_of)
END IF
END SUBROUTINE write_forces
END MODULE qs_force