-
Notifications
You must be signed in to change notification settings - Fork 1
/
mixed_cdft_methods.F
4210 lines (4101 loc) · 237 KB
/
mixed_cdft_methods.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
684
685
686
687
688
689
690
691
692
693
694
695
696
697
698
699
700
701
702
703
704
705
706
707
708
709
710
711
712
713
714
715
716
717
718
719
720
721
722
723
724
725
726
727
728
729
730
731
732
733
734
735
736
737
738
739
740
741
742
743
744
745
746
747
748
749
750
751
752
753
754
755
756
757
758
759
760
761
762
763
764
765
766
767
768
769
770
771
772
773
774
775
776
777
778
779
780
781
782
783
784
785
786
787
788
789
790
791
792
793
794
795
796
797
798
799
800
801
802
803
804
805
806
807
808
809
810
811
812
813
814
815
816
817
818
819
820
821
822
823
824
825
826
827
828
829
830
831
832
833
834
835
836
837
838
839
840
841
842
843
844
845
846
847
848
849
850
851
852
853
854
855
856
857
858
859
860
861
862
863
864
865
866
867
868
869
870
871
872
873
874
875
876
877
878
879
880
881
882
883
884
885
886
887
888
889
890
891
892
893
894
895
896
897
898
899
900
901
902
903
904
905
906
907
908
909
910
911
912
913
914
915
916
917
918
919
920
921
922
923
924
925
926
927
928
929
930
931
932
933
934
935
936
937
938
939
940
941
942
943
944
945
946
947
948
949
950
951
952
953
954
955
956
957
958
959
960
961
962
963
964
965
966
967
968
969
970
971
972
973
974
975
976
977
978
979
980
981
982
983
984
985
986
987
988
989
990
991
992
993
994
995
996
997
998
999
1000
!--------------------------------------------------------------------------------------------------!
! 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 Methods for mixed CDFT calculations
!> \par History
!> Separated CDFT routines from mixed_environment_utils
!> \author Nico Holmberg [01.2017]
! **************************************************************************************************
MODULE mixed_cdft_methods
USE ao_util, ONLY: exp_radius_very_extended
USE atomic_kind_types, ONLY: atomic_kind_type,&
get_atomic_kind
USE cell_types, ONLY: cell_type,&
pbc
USE cp_array_utils, ONLY: cp_1d_i_p_type,&
cp_1d_r_p_type,&
cp_2d_r_p_type
USE cp_control_types, ONLY: dft_control_type
USE cp_dbcsr_api, ONLY: &
dbcsr_add, dbcsr_copy, dbcsr_create, dbcsr_init_p, dbcsr_p_type, dbcsr_release, &
dbcsr_release_p, dbcsr_scale, dbcsr_type
USE cp_dbcsr_diag, ONLY: cp_dbcsr_syevd
USE cp_dbcsr_operations, ONLY: cp_dbcsr_sm_fm_multiply
USE cp_fm_basic_linalg, ONLY: cp_fm_column_scale,&
cp_fm_invert,&
cp_fm_transpose
USE cp_fm_struct, ONLY: cp_fm_struct_create,&
cp_fm_struct_release,&
cp_fm_struct_type
USE cp_fm_types, ONLY: cp_fm_create,&
cp_fm_get_info,&
cp_fm_release,&
cp_fm_set_all,&
cp_fm_to_fm,&
cp_fm_type,&
cp_fm_write_formatted
USE cp_log_handling, ONLY: cp_get_default_logger,&
cp_logger_type,&
cp_to_string
USE cp_output_handling, ONLY: cp_print_key_finished_output,&
cp_print_key_unit_nr
USE cp_realspace_grid_cube, ONLY: cp_pw_to_cube
USE cp_subsys_types, ONLY: cp_subsys_get,&
cp_subsys_type
USE cp_units, ONLY: cp_unit_from_cp2k
USE force_env_types, ONLY: force_env_get,&
force_env_type,&
use_qmmm,&
use_qmmmx,&
use_qs_force
USE grid_api, ONLY: GRID_FUNC_AB,&
collocate_pgf_product
USE hirshfeld_methods, ONLY: create_shape_function
USE hirshfeld_types, ONLY: hirshfeld_type
USE input_constants, ONLY: &
becke_cutoff_element, becke_cutoff_global, cdft_alpha_constraint, cdft_beta_constraint, &
cdft_charge_constraint, cdft_magnetization_constraint, mix_cdft, mixed_cdft_parallel, &
mixed_cdft_parallel_nobuild, mixed_cdft_serial, outer_scf_becke_constraint
USE input_section_types, ONLY: section_get_lval,&
section_vals_get,&
section_vals_get_subs_vals,&
section_vals_type,&
section_vals_val_get
USE kinds, ONLY: default_path_length,&
dp
USE machine, ONLY: m_walltime
USE mathlib, ONLY: diamat_all
USE memory_utilities, ONLY: reallocate
USE message_passing, ONLY: mp_request_type,&
mp_testall,&
mp_waitall
USE mixed_cdft_types, ONLY: mixed_cdft_result_type_set,&
mixed_cdft_settings_type,&
mixed_cdft_type,&
mixed_cdft_type_create,&
mixed_cdft_work_type_release
USE mixed_cdft_utils, ONLY: &
hfun_zero, map_permutation_to_states, mixed_cdft_assemble_block_diag, &
mixed_cdft_diagonalize_blocks, mixed_cdft_get_blocks, mixed_cdft_init_structures, &
mixed_cdft_parse_settings, mixed_cdft_print_couplings, mixed_cdft_read_block_diag, &
mixed_cdft_redistribute_arrays, mixed_cdft_release_work, mixed_cdft_transfer_settings
USE mixed_environment_types, ONLY: get_mixed_env,&
mixed_environment_type,&
set_mixed_env
USE parallel_gemm_api, ONLY: parallel_gemm
USE particle_list_types, ONLY: particle_list_type
USE particle_types, ONLY: particle_type
USE pw_env_types, ONLY: pw_env_get,&
pw_env_type
USE pw_methods, ONLY: pw_copy,&
pw_scale,&
pw_zero
USE pw_pool_types, ONLY: pw_pool_type
USE qs_cdft_types, ONLY: cdft_control_type
USE qs_energy_types, ONLY: qs_energy_type
USE qs_environment_types, ONLY: get_qs_env,&
qs_environment_type
USE qs_kind_types, ONLY: qs_kind_type
USE qs_mo_io, ONLY: read_mo_set_from_restart,&
wfn_restart_file_name
USE qs_mo_methods, ONLY: make_basis_simple,&
make_basis_sm
USE qs_mo_types, ONLY: allocate_mo_set,&
deallocate_mo_set,&
mo_set_type,&
set_mo_set
USE realspace_grid_types, ONLY: realspace_grid_type,&
rs_grid_zero,&
transfer_rs2pw
USE util, ONLY: sort
#include "./base/base_uses.f90"
IMPLICIT NONE
PRIVATE
CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'mixed_cdft_methods'
LOGICAL, PARAMETER, PRIVATE :: debug_this_module = .FALSE.
PUBLIC :: mixed_cdft_init, &
mixed_cdft_build_weight, &
mixed_cdft_calculate_coupling
CONTAINS
! **************************************************************************************************
!> \brief Initialize a mixed CDFT calculation
!> \param force_env the force_env that holds the CDFT states
!> \param calculate_forces determines if forces should be calculated
!> \par History
!> 01.2016 created [Nico Holmberg]
! **************************************************************************************************
SUBROUTINE mixed_cdft_init(force_env, calculate_forces)
TYPE(force_env_type), POINTER :: force_env
LOGICAL, INTENT(IN) :: calculate_forces
CHARACTER(len=*), PARAMETER :: routineN = 'mixed_cdft_init'
INTEGER :: et_freq, handle, iforce_eval, iounit, &
mixing_type, nforce_eval
LOGICAL :: explicit, is_parallel, is_qmmm
TYPE(cp_logger_type), POINTER :: logger
TYPE(cp_subsys_type), POINTER :: subsys_mix
TYPE(force_env_type), POINTER :: force_env_qs
TYPE(mixed_cdft_settings_type) :: settings
TYPE(mixed_cdft_type), POINTER :: mixed_cdft
TYPE(mixed_environment_type), POINTER :: mixed_env
TYPE(particle_list_type), POINTER :: particles_mix
TYPE(section_vals_type), POINTER :: force_env_section, mapping_section, &
md_section, mixed_section, &
print_section, root_section
NULLIFY (subsys_mix, force_env_qs, force_env_section, print_section, &
root_section, mixed_section, md_section, mixed_env, mixed_cdft, &
mapping_section)
NULLIFY (settings%grid_span, settings%npts, settings%cutoff, settings%rel_cutoff, &
settings%spherical, settings%rs_dims, settings%odd, settings%atoms, &
settings%coeffs, settings%si, settings%sr, &
settings%cutoffs, settings%radii)
is_qmmm = .FALSE.
logger => cp_get_default_logger()
CPASSERT(ASSOCIATED(force_env))
nforce_eval = SIZE(force_env%sub_force_env)
CALL timeset(routineN, handle)
CALL force_env_get(force_env=force_env, force_env_section=force_env_section)
mixed_env => force_env%mixed_env
mixed_section => section_vals_get_subs_vals(force_env_section, "MIXED")
print_section => section_vals_get_subs_vals(force_env_section, "MIXED%MIXED_CDFT%PRINT%PROGRAM_RUN_INFO")
iounit = cp_print_key_unit_nr(logger, print_section, '', extension='.mixedLog')
! Check if a mixed CDFT calculation is requested
CALL section_vals_val_get(mixed_section, "MIXING_TYPE", i_val=mixing_type)
IF (mixing_type == mix_cdft .AND. .NOT. ASSOCIATED(mixed_env%cdft_control)) THEN
mixed_env%do_mixed_cdft = .TRUE.
IF (mixed_env%do_mixed_cdft) THEN
! Sanity check
IF (nforce_eval .LT. 2) &
CALL cp_abort(__LOCATION__, &
"Mixed CDFT calculation requires at least 2 force_evals.")
mapping_section => section_vals_get_subs_vals(mixed_section, "MAPPING")
CALL section_vals_get(mapping_section, explicit=explicit)
! The sub_force_envs must share the same geometrical structure
IF (explicit) &
CPABORT("Please disable section &MAPPING for mixed CDFT calculations")
CALL section_vals_val_get(mixed_section, "MIXED_CDFT%COUPLING", i_val=et_freq)
IF (et_freq .LT. 0) THEN
mixed_env%do_mixed_et = .FALSE.
ELSE
mixed_env%do_mixed_et = .TRUE.
IF (et_freq == 0) THEN
mixed_env%et_freq = 1
ELSE
mixed_env%et_freq = et_freq
END IF
END IF
! Start initializing the mixed_cdft type
! First determine if the calculation is pure DFT or QMMM and find the qs force_env
DO iforce_eval = 1, nforce_eval
IF (.NOT. ASSOCIATED(force_env%sub_force_env(iforce_eval)%force_env)) CYCLE
SELECT CASE (force_env%sub_force_env(iforce_eval)%force_env%in_use)
CASE (use_qs_force)
force_env_qs => force_env%sub_force_env(iforce_eval)%force_env
CASE (use_qmmm)
is_qmmm = .TRUE.
! This is really the container for QMMM
force_env_qs => force_env%sub_force_env(iforce_eval)%force_env
CASE (use_qmmmx)
CPABORT("No force mixing allowed for mixed CDFT QM/MM")
CASE DEFAULT
CPASSERT(.FALSE.)
END SELECT
CPASSERT(ASSOCIATED(force_env_qs))
END DO
! Get infos about the mixed subsys
IF (.NOT. is_qmmm) THEN
CALL force_env_get(force_env=force_env, &
subsys=subsys_mix)
CALL cp_subsys_get(subsys=subsys_mix, &
particles=particles_mix)
ELSE
CALL get_qs_env(force_env_qs%qmmm_env%qs_env, &
cp_subsys=subsys_mix)
CALL cp_subsys_get(subsys=subsys_mix, &
particles=particles_mix)
END IF
! Init mixed_cdft_type
ALLOCATE (mixed_cdft)
CALL mixed_cdft_type_create(mixed_cdft)
mixed_cdft%first_iteration = .TRUE.
! Determine what run type to use
IF (mixed_env%ngroups == 1) THEN
! States treated in serial, possibly copying CDFT weight function and gradients from state to state
mixed_cdft%run_type = mixed_cdft_serial
ELSE IF (mixed_env%ngroups == 2) THEN
CALL section_vals_val_get(mixed_section, "MIXED_CDFT%PARALLEL_BUILD", l_val=is_parallel)
IF (is_parallel) THEN
! Treat states in parallel, build weight function and gradients in parallel before SCF process
mixed_cdft%run_type = mixed_cdft_parallel
IF (.NOT. nforce_eval == 2) &
CALL cp_abort(__LOCATION__, &
"Parallel mode mixed CDFT calculation supports only 2 force_evals.")
ELSE
! Treat states in parallel, but each states builds its own weight function and gradients
mixed_cdft%run_type = mixed_cdft_parallel_nobuild
END IF
ELSE
mixed_cdft%run_type = mixed_cdft_parallel_nobuild
END IF
! Store QMMM flag
mixed_env%do_mixed_qmmm_cdft = is_qmmm
! Setup dynamic load balancing
CALL section_vals_val_get(mixed_section, "MIXED_CDFT%DLB", l_val=mixed_cdft%dlb)
mixed_cdft%dlb = mixed_cdft%dlb .AND. calculate_forces ! disable if forces are not needed
mixed_cdft%dlb = mixed_cdft%dlb .AND. (mixed_cdft%run_type == mixed_cdft_parallel) ! disable if not parallel
IF (mixed_cdft%dlb) THEN
ALLOCATE (mixed_cdft%dlb_control)
NULLIFY (mixed_cdft%dlb_control%weight, mixed_cdft%dlb_control%gradients, &
mixed_cdft%dlb_control%cavity, mixed_cdft%dlb_control%target_list, &
mixed_cdft%dlb_control%bo, mixed_cdft%dlb_control%expected_work, &
mixed_cdft%dlb_control%prediction_error, mixed_cdft%dlb_control%sendbuff, &
mixed_cdft%dlb_control%recvbuff, mixed_cdft%dlb_control%recv_work_repl, &
mixed_cdft%dlb_control%recv_info)
CALL section_vals_val_get(mixed_section, "MIXED_CDFT%LOAD_SCALE", &
r_val=mixed_cdft%dlb_control%load_scale)
CALL section_vals_val_get(mixed_section, "MIXED_CDFT%VERY_OVERLOADED", &
r_val=mixed_cdft%dlb_control%very_overloaded)
CALL section_vals_val_get(mixed_section, "MIXED_CDFT%MORE_WORK", &
i_val=mixed_cdft%dlb_control%more_work)
END IF
! Metric/Wavefunction overlap method/Lowdin orthogonalization/CDFT-CI
mixed_cdft%calculate_metric = .FALSE.
mixed_cdft%wfn_overlap_method = .FALSE.
mixed_cdft%use_lowdin = .FALSE.
mixed_cdft%do_ci = .FALSE.
mixed_cdft%nonortho_coupling = .FALSE.
mixed_cdft%identical_constraints = .TRUE.
mixed_cdft%block_diagonalize = .FALSE.
IF (mixed_env%do_mixed_et) THEN
CALL section_vals_val_get(mixed_section, "MIXED_CDFT%METRIC", &
l_val=mixed_cdft%calculate_metric)
CALL section_vals_val_get(mixed_section, "MIXED_CDFT%WFN_OVERLAP", &
l_val=mixed_cdft%wfn_overlap_method)
CALL section_vals_val_get(mixed_section, "MIXED_CDFT%LOWDIN", &
l_val=mixed_cdft%use_lowdin)
CALL section_vals_val_get(mixed_section, "MIXED_CDFT%CI", &
l_val=mixed_cdft%do_ci)
CALL section_vals_val_get(mixed_section, "MIXED_CDFT%NONORTHOGONAL_COUPLING", &
l_val=mixed_cdft%nonortho_coupling)
CALL section_vals_val_get(mixed_section, "MIXED_CDFT%BLOCK_DIAGONALIZE", &
l_val=mixed_cdft%block_diagonalize)
END IF
! Inversion method
CALL section_vals_val_get(mixed_section, "MIXED_CDFT%EPS_SVD", r_val=mixed_cdft%eps_svd)
IF (mixed_cdft%eps_svd .LT. 0.0_dp .OR. mixed_cdft%eps_svd .GT. 1.0_dp) &
CPABORT("Illegal value for EPS_SVD. Value must be between 0.0 and 1.0.")
! MD related settings
CALL force_env_get(force_env, root_section=root_section)
md_section => section_vals_get_subs_vals(root_section, "MOTION%MD")
CALL section_vals_val_get(md_section, "TIMESTEP", r_val=mixed_cdft%sim_dt)
CALL section_vals_val_get(md_section, "STEP_START_VAL", i_val=mixed_cdft%sim_step)
mixed_cdft%sim_step = mixed_cdft%sim_step - 1 ! to get the first step correct
mixed_cdft%sim_dt = cp_unit_from_cp2k(mixed_cdft%sim_dt, "fs")
! Parse constraint settings from the individual force_evals and check consistency
CALL mixed_cdft_parse_settings(force_env, mixed_env, mixed_cdft, &
settings, natom=SIZE(particles_mix%els))
! Transfer settings to mixed_cdft
CALL mixed_cdft_transfer_settings(force_env, mixed_cdft, settings)
! Initilize necessary structures
CALL mixed_cdft_init_structures(force_env, force_env_qs, mixed_env, mixed_cdft, settings)
! Write information about the mixed CDFT calculation
IF (iounit > 0) THEN
WRITE (iounit, *) ""
WRITE (iounit, FMT="(T2,A,T71)") &
"MIXED_CDFT| Activating mixed CDFT calculation"
WRITE (iounit, FMT="(T2,A,T71,I10)") &
"MIXED_CDFT| Number of CDFT states: ", nforce_eval
SELECT CASE (mixed_cdft%run_type)
CASE (mixed_cdft_parallel)
WRITE (iounit, FMT="(T2,A,T71)") &
"MIXED_CDFT| CDFT states calculation mode: parallel with build"
WRITE (iounit, FMT="(T2,A,T71)") &
"MIXED_CDFT| Becke constraint is first built using all available processors"
WRITE (iounit, FMT="(T2,A,T71)") &
" and then copied to both states with their own processor groups"
CASE (mixed_cdft_serial)
WRITE (iounit, FMT="(T2,A,T71)") &
"MIXED_CDFT| CDFT states calculation mode: serial"
IF (mixed_cdft%identical_constraints) THEN
WRITE (iounit, FMT="(T2,A,T71)") &
"MIXED_CDFT| The constraints are built before the SCF procedure of the first"
WRITE (iounit, FMT="(T2,A,T71)") &
" CDFT state and subsequently copied to the other states"
ELSE
WRITE (iounit, FMT="(T2,A,T71)") &
"MIXED_CDFT| The constraints are separately built for all CDFT states"
END IF
CASE (mixed_cdft_parallel_nobuild)
WRITE (iounit, FMT="(T2,A,T71)") &
"MIXED_CDFT| CDFT states calculation mode: parallel without build"
WRITE (iounit, FMT="(T2,A,T71)") &
"MIXED_CDFT| The constraints are separately built for all CDFT states"
CASE DEFAULT
CPABORT("Unknown mixed CDFT run type.")
END SELECT
WRITE (iounit, FMT="(T2,A,T71,L10)") &
"MIXED_CDFT| Calculating electronic coupling between states: ", mixed_env%do_mixed_et
WRITE (iounit, FMT="(T2,A,T71,L10)") &
"MIXED_CDFT| Calculating electronic coupling reliability metric: ", mixed_cdft%calculate_metric
WRITE (iounit, FMT="(T2,A,T71,L10)") &
"MIXED_CDFT| Configuration interaction (CDFT-CI) was requested: ", mixed_cdft%do_ci
WRITE (iounit, FMT="(T2,A,T71,L10)") &
"MIXED_CDFT| Block diagonalizing the mixed CDFT Hamiltonian: ", mixed_cdft%block_diagonalize
IF (mixed_cdft%run_type == mixed_cdft_parallel) THEN
WRITE (iounit, FMT="(T2,A,T71,L10)") &
"MIXED_CDFT| Dynamic load balancing enabled: ", mixed_cdft%dlb
IF (mixed_cdft%dlb) THEN
WRITE (iounit, FMT="(T2,A,T71)") "MIXED_CDFT| Dynamic load balancing parameters:"
WRITE (iounit, FMT="(T2,A,T71,F10.2)") &
"MIXED_CDFT| load_scale:", mixed_cdft%dlb_control%load_scale
WRITE (iounit, FMT="(T2,A,T71,F10.2)") &
"MIXED_CDFT| very_overloaded:", mixed_cdft%dlb_control%very_overloaded
WRITE (iounit, FMT="(T2,A,T71,I10)") &
"MIXED_CDFT| more_work:", mixed_cdft%dlb_control%more_work
END IF
END IF
IF (mixed_env%do_mixed_et) THEN
IF (mixed_cdft%eps_svd == 0.0_dp) THEN
WRITE (iounit, FMT="(T2,A,T71)") "MIXED_CDFT| Matrix inversions calculated with LU decomposition."
ELSE
WRITE (iounit, FMT="(T2,A,T71)") "MIXED_CDFT| Matrix inversions calculated with SVD decomposition."
WRITE (iounit, FMT="(T2,A,T71,ES10.2)") "MIXED_CDFT| EPS_SVD:", mixed_cdft%eps_svd
END IF
END IF
END IF
CALL set_mixed_env(mixed_env, cdft_control=mixed_cdft)
END IF
END IF
CALL cp_print_key_finished_output(iounit, logger, force_env_section, &
"MIXED%MIXED_CDFT%PRINT%PROGRAM_RUN_INFO")
CALL timestop(handle)
END SUBROUTINE mixed_cdft_init
! **************************************************************************************************
!> \brief Driver routine to handle the build of CDFT weight/gradient in parallel and serial modes
!> \param force_env the force_env that holds the CDFT states
!> \param calculate_forces if forces should be calculated
!> \param iforce_eval index of the currently active CDFT state (serial mode only)
!> \par History
!> 01.2017 created [Nico Holmberg]
! **************************************************************************************************
SUBROUTINE mixed_cdft_build_weight(force_env, calculate_forces, iforce_eval)
TYPE(force_env_type), POINTER :: force_env
LOGICAL, INTENT(IN) :: calculate_forces
INTEGER, INTENT(IN), OPTIONAL :: iforce_eval
TYPE(mixed_cdft_type), POINTER :: mixed_cdft
NULLIFY (mixed_cdft)
CPASSERT(ASSOCIATED(force_env))
CALL get_mixed_env(force_env%mixed_env, cdft_control=mixed_cdft)
CPASSERT(ASSOCIATED(mixed_cdft))
IF (.NOT. PRESENT(iforce_eval)) THEN
SELECT CASE (mixed_cdft%run_type)
CASE (mixed_cdft_parallel)
CALL mixed_cdft_build_weight_parallel(force_env, calculate_forces)
CASE (mixed_cdft_parallel_nobuild)
CALL mixed_cdft_set_flags(force_env)
CASE DEFAULT
! Do nothing
END SELECT
ELSE
SELECT CASE (mixed_cdft%run_type)
CASE (mixed_cdft_serial)
CALL mixed_cdft_transfer_weight(force_env, calculate_forces, iforce_eval)
CASE DEFAULT
! Do nothing
END SELECT
END IF
END SUBROUTINE mixed_cdft_build_weight
! **************************************************************************************************
!> \brief Build CDFT weight and gradient on 2N processors and copy it to two N processor subgroups
!> \param force_env the force_env that holds the CDFT states
!> \param calculate_forces if forces should be calculted
!> \par History
!> 01.2016 created [Nico Holmberg]
! **************************************************************************************************
SUBROUTINE mixed_cdft_build_weight_parallel(force_env, calculate_forces)
TYPE(force_env_type), POINTER :: force_env
LOGICAL, INTENT(IN) :: calculate_forces
CHARACTER(len=*), PARAMETER :: routineN = 'mixed_cdft_build_weight_parallel'
INTEGER :: handle, handle2, i, iforce_eval, ind, INDEX(6), iounit, j, lb_min, &
my_special_work, natom, nforce_eval, recv_offset, ub_max
INTEGER, DIMENSION(2, 3) :: bo
INTEGER, DIMENSION(:), POINTER :: lb, sendbuffer_i, ub
REAL(KIND=dp) :: t1, t2
TYPE(cdft_control_type), POINTER :: cdft_control, cdft_control_target
TYPE(cp_logger_type), POINTER :: logger
TYPE(cp_subsys_type), POINTER :: subsys_mix
TYPE(dft_control_type), POINTER :: dft_control
TYPE(force_env_type), POINTER :: force_env_qs
TYPE(mixed_cdft_type), POINTER :: mixed_cdft
TYPE(mixed_environment_type), POINTER :: mixed_env
TYPE(mp_request_type), DIMENSION(:), POINTER :: req_total
TYPE(particle_list_type), POINTER :: particles_mix
TYPE(pw_env_type), POINTER :: pw_env
TYPE(pw_pool_type), POINTER :: auxbas_pw_pool, mixed_auxbas_pw_pool
TYPE(qs_environment_type), POINTER :: qs_env
TYPE(section_vals_type), POINTER :: force_env_section, print_section
TYPE buffers
INTEGER :: imap(6)
INTEGER, DIMENSION(:), &
POINTER :: iv => null()
REAL(KIND=dp), POINTER, &
DIMENSION(:, :, :) :: r3 => null()
REAL(KIND=dp), POINTER, &
DIMENSION(:, :, :, :) :: r4 => null()
END TYPE buffers
TYPE(buffers), DIMENSION(:), POINTER :: recvbuffer
NULLIFY (subsys_mix, force_env_qs, particles_mix, force_env_section, print_section, &
mixed_env, mixed_cdft, pw_env, auxbas_pw_pool, mixed_auxbas_pw_pool, &
qs_env, dft_control, sendbuffer_i, lb, ub, req_total, recvbuffer, &
cdft_control, cdft_control_target)
logger => cp_get_default_logger()
CPASSERT(ASSOCIATED(force_env))
nforce_eval = SIZE(force_env%sub_force_env)
CALL timeset(routineN, handle)
t1 = m_walltime()
! Get infos about the mixed subsys
CALL force_env_get(force_env=force_env, &
subsys=subsys_mix, &
force_env_section=force_env_section)
CALL cp_subsys_get(subsys=subsys_mix, &
particles=particles_mix)
DO iforce_eval = 1, nforce_eval
IF (.NOT. ASSOCIATED(force_env%sub_force_env(iforce_eval)%force_env)) CYCLE
SELECT CASE (force_env%sub_force_env(iforce_eval)%force_env%in_use)
CASE (use_qs_force)
force_env_qs => force_env%sub_force_env(iforce_eval)%force_env
CASE (use_qmmm)
force_env_qs => force_env%sub_force_env(iforce_eval)%force_env
CASE DEFAULT
CPASSERT(.FALSE.)
END SELECT
END DO
IF (.NOT. force_env%mixed_env%do_mixed_qmmm_cdft) THEN
CALL force_env_get(force_env=force_env_qs, &
qs_env=qs_env, &
subsys=subsys_mix)
CALL cp_subsys_get(subsys=subsys_mix, &
particles=particles_mix)
ELSE
qs_env => force_env_qs%qmmm_env%qs_env
CALL get_qs_env(qs_env, cp_subsys=subsys_mix)
CALL cp_subsys_get(subsys=subsys_mix, &
particles=particles_mix)
END IF
mixed_env => force_env%mixed_env
print_section => section_vals_get_subs_vals(force_env_section, "MIXED%MIXED_CDFT%PRINT%PROGRAM_RUN_INFO")
iounit = cp_print_key_unit_nr(logger, print_section, '', extension='.mixedLog')
CALL get_mixed_env(mixed_env, cdft_control=mixed_cdft)
CPASSERT(ASSOCIATED(mixed_cdft))
cdft_control => mixed_cdft%cdft_control
CPASSERT(ASSOCIATED(cdft_control))
! Calculate the Becke weight function and gradient on all np processors
CALL pw_env_get(pw_env=mixed_cdft%pw_env, auxbas_pw_pool=mixed_auxbas_pw_pool)
natom = SIZE(particles_mix%els)
CALL mixed_becke_constraint(force_env, calculate_forces)
! Start replicating the working arrays on both np/2 processor groups
mixed_cdft%sim_step = mixed_cdft%sim_step + 1
CALL get_qs_env(qs_env, pw_env=pw_env, dft_control=dft_control)
cdft_control_target => dft_control%qs_control%cdft_control
CPASSERT(dft_control%qs_control%cdft)
CPASSERT(ASSOCIATED(cdft_control_target))
CALL pw_env_get(pw_env, auxbas_pw_pool=auxbas_pw_pool)
bo = auxbas_pw_pool%pw_grid%bounds_local
! First determine the size of the arrays along the confinement dir
IF (mixed_cdft%is_special) THEN
my_special_work = 2
ELSE
my_special_work = 1
END IF
ALLOCATE (recvbuffer(SIZE(mixed_cdft%source_list)))
ALLOCATE (req_total(my_special_work*SIZE(mixed_cdft%dest_list) + SIZE(mixed_cdft%source_list)))
ALLOCATE (lb(SIZE(mixed_cdft%source_list)), ub(SIZE(mixed_cdft%source_list)))
IF (cdft_control%becke_control%cavity_confine) THEN
! Gaussian confinement => the bounds depend on the processor and need to be communicated
ALLOCATE (sendbuffer_i(2))
sendbuffer_i = cdft_control%becke_control%confine_bounds
DO i = 1, SIZE(mixed_cdft%source_list)
ALLOCATE (recvbuffer(i)%iv(2))
CALL force_env%para_env%irecv(msgout=recvbuffer(i)%iv, source=mixed_cdft%source_list(i), &
request=req_total(i))
END DO
DO i = 1, my_special_work
DO j = 1, SIZE(mixed_cdft%dest_list)
ind = j + (i - 1)*SIZE(mixed_cdft%dest_list) + SIZE(mixed_cdft%source_list)
CALL force_env%para_env%isend(msgin=sendbuffer_i, &
dest=mixed_cdft%dest_list(j) + (i - 1)*force_env%para_env%num_pe/2, &
request=req_total(ind))
END DO
END DO
CALL mp_waitall(req_total)
! Find the largest/smallest value of ub/lb
DEALLOCATE (sendbuffer_i)
lb_min = HUGE(0)
ub_max = -HUGE(0)
DO i = 1, SIZE(mixed_cdft%source_list)
lb(i) = recvbuffer(i)%iv(1)
ub(i) = recvbuffer(i)%iv(2)
IF (lb(i) .LT. lb_min) lb_min = lb(i)
IF (ub(i) .GT. ub_max) ub_max = ub(i)
DEALLOCATE (recvbuffer(i)%iv)
END DO
! Take into account the grids already communicated during dlb
IF (mixed_cdft%dlb) THEN
IF (ANY(mixed_cdft%dlb_control%recv_work_repl)) THEN
DO j = 1, SIZE(mixed_cdft%dlb_control%recv_work_repl)
IF (mixed_cdft%dlb_control%recv_work_repl(j)) THEN
DO i = 1, SIZE(mixed_cdft%dlb_control%recvbuff(j)%buffs)
IF (LBOUND(mixed_cdft%dlb_control%recvbuff(j)%buffs(i)%weight, 3) &
.LT. lb_min) lb_min = LBOUND(mixed_cdft%dlb_control%recvbuff(j)%buffs(i)%weight, 3)
IF (UBOUND(mixed_cdft%dlb_control%recvbuff(j)%buffs(i)%weight, 3) &
.GT. ub_max) ub_max = UBOUND(mixed_cdft%dlb_control%recvbuff(j)%buffs(i)%weight, 3)
END DO
END IF
END DO
END IF
END IF
ELSE
! No confinement
ub_max = bo(2, 3)
lb_min = bo(1, 3)
lb = lb_min
ub = ub_max
END IF
! Determine the sender specific indices of grid slices that are to be received
CALL timeset(routineN//"_comm", handle2)
DO j = 1, SIZE(recvbuffer)
ind = j + (j/2)
IF (mixed_cdft%is_special) THEN
recvbuffer(j)%imap = (/mixed_cdft%source_list_bo(1, j), mixed_cdft%source_list_bo(2, j), &
mixed_cdft%source_list_bo(3, j), mixed_cdft%source_list_bo(4, j), &
lb(j), ub(j)/)
ELSE IF (mixed_cdft%is_pencil) THEN
recvbuffer(j)%imap = (/bo(1, 1), bo(2, 1), mixed_cdft%recv_bo(ind), mixed_cdft%recv_bo(ind + 1), lb(j), ub(j)/)
ELSE
recvbuffer(j)%imap = (/mixed_cdft%recv_bo(ind), mixed_cdft%recv_bo(ind + 1), bo(1, 2), bo(2, 2), lb(j), ub(j)/)
END IF
END DO
IF (mixed_cdft%dlb .AND. .NOT. mixed_cdft%is_special) THEN
IF (mixed_cdft%dlb_control%recv_work_repl(1) .OR. mixed_cdft%dlb_control%recv_work_repl(2)) THEN
DO j = 1, 2
recv_offset = 0
IF (mixed_cdft%dlb_control%recv_work_repl(j)) &
recv_offset = SUM(mixed_cdft%dlb_control%recv_info(j)%target_list(2, :))
IF (mixed_cdft%is_pencil) THEN
recvbuffer(j)%imap(1) = recvbuffer(j)%imap(1) + recv_offset
ELSE
recvbuffer(j)%imap(3) = recvbuffer(j)%imap(3) + recv_offset
END IF
END DO
END IF
END IF
! Transfer the arrays one-by-one and deallocate shared storage
! Start with the weight function
DO j = 1, SIZE(mixed_cdft%source_list)
ALLOCATE (recvbuffer(j)%r3(recvbuffer(j)%imap(1):recvbuffer(j)%imap(2), &
recvbuffer(j)%imap(3):recvbuffer(j)%imap(4), &
recvbuffer(j)%imap(5):recvbuffer(j)%imap(6)))
CALL force_env%para_env%irecv(msgout=recvbuffer(j)%r3, source=mixed_cdft%source_list(j), &
request=req_total(j))
END DO
DO i = 1, my_special_work
DO j = 1, SIZE(mixed_cdft%dest_list)
ind = j + (i - 1)*SIZE(mixed_cdft%dest_list) + SIZE(mixed_cdft%source_list)
IF (mixed_cdft%is_special) THEN
CALL force_env%para_env%isend(msgin=mixed_cdft%sendbuff(j)%weight, &
dest=mixed_cdft%dest_list(j) + (i - 1)*force_env%para_env%num_pe/2, &
request=req_total(ind))
ELSE
CALL force_env%para_env%isend(msgin=mixed_cdft%weight, dest=mixed_cdft%dest_list(j), &
request=req_total(ind))
END IF
END DO
END DO
CALL mp_waitall(req_total)
IF (mixed_cdft%is_special) THEN
DO j = 1, SIZE(mixed_cdft%dest_list)
DEALLOCATE (mixed_cdft%sendbuff(j)%weight)
END DO
ELSE
DEALLOCATE (mixed_cdft%weight)
END IF
! In principle, we could reduce the memory footprint of becke_pot by only transferring the nonzero portion
! of the potential, but this would require a custom integrate_v_rspace
ALLOCATE (cdft_control_target%group(1)%weight)
CALL auxbas_pw_pool%create_pw(cdft_control_target%group(1)%weight)
CALL pw_zero(cdft_control_target%group(1)%weight)
! Assemble the recved slices
DO j = 1, SIZE(mixed_cdft%source_list)
cdft_control_target%group(1)%weight%array(recvbuffer(j)%imap(1):recvbuffer(j)%imap(2), &
recvbuffer(j)%imap(3):recvbuffer(j)%imap(4), &
recvbuffer(j)%imap(5):recvbuffer(j)%imap(6)) = recvbuffer(j)%r3
END DO
! Do the same for slices sent during dlb
IF (mixed_cdft%dlb) THEN
IF (ANY(mixed_cdft%dlb_control%recv_work_repl)) THEN
DO j = 1, SIZE(mixed_cdft%dlb_control%recv_work_repl)
IF (mixed_cdft%dlb_control%recv_work_repl(j)) THEN
DO i = 1, SIZE(mixed_cdft%dlb_control%recvbuff(j)%buffs)
index = (/LBOUND(mixed_cdft%dlb_control%recvbuff(j)%buffs(i)%weight, 1), &
UBOUND(mixed_cdft%dlb_control%recvbuff(j)%buffs(i)%weight, 1), &
LBOUND(mixed_cdft%dlb_control%recvbuff(j)%buffs(i)%weight, 2), &
UBOUND(mixed_cdft%dlb_control%recvbuff(j)%buffs(i)%weight, 2), &
LBOUND(mixed_cdft%dlb_control%recvbuff(j)%buffs(i)%weight, 3), &
UBOUND(mixed_cdft%dlb_control%recvbuff(j)%buffs(i)%weight, 3)/)
cdft_control_target%group(1)%weight%array(INDEX(1):INDEX(2), &
INDEX(3):INDEX(4), &
INDEX(5):INDEX(6)) = &
mixed_cdft%dlb_control%recvbuff(j)%buffs(i)%weight
DEALLOCATE (mixed_cdft%dlb_control%recvbuff(j)%buffs(i)%weight)
END DO
END IF
END DO
END IF
END IF
! Gaussian confinement cavity
IF (cdft_control%becke_control%cavity_confine) THEN
DO j = 1, SIZE(mixed_cdft%source_list)
CALL force_env%para_env%irecv(msgout=recvbuffer(j)%r3, source=mixed_cdft%source_list(j), &
request=req_total(j))
END DO
DO i = 1, my_special_work
DO j = 1, SIZE(mixed_cdft%dest_list)
ind = j + (i - 1)*SIZE(mixed_cdft%dest_list) + SIZE(mixed_cdft%source_list)
IF (mixed_cdft%is_special) THEN
CALL force_env%para_env%isend(msgin=mixed_cdft%sendbuff(j)%cavity, &
dest=mixed_cdft%dest_list(j) + (i - 1)*force_env%para_env%num_pe/2, &
request=req_total(ind))
ELSE
CALL force_env%para_env%isend(msgin=mixed_cdft%cavity, dest=mixed_cdft%dest_list(j), &
request=req_total(ind))
END IF
END DO
END DO
CALL mp_waitall(req_total)
IF (mixed_cdft%is_special) THEN
DO j = 1, SIZE(mixed_cdft%dest_list)
DEALLOCATE (mixed_cdft%sendbuff(j)%cavity)
END DO
ELSE
DEALLOCATE (mixed_cdft%cavity)
END IF
! We only need the nonzero part of the confinement cavity
ALLOCATE (cdft_control_target%becke_control%cavity_mat(bo(1, 1):bo(2, 1), &
bo(1, 2):bo(2, 2), &
lb_min:ub_max))
cdft_control_target%becke_control%cavity_mat = 0.0_dp
DO j = 1, SIZE(mixed_cdft%source_list)
cdft_control_target%becke_control%cavity_mat(recvbuffer(j)%imap(1):recvbuffer(j)%imap(2), &
recvbuffer(j)%imap(3):recvbuffer(j)%imap(4), &
recvbuffer(j)%imap(5):recvbuffer(j)%imap(6)) = recvbuffer(j)%r3
END DO
IF (mixed_cdft%dlb) THEN
IF (ANY(mixed_cdft%dlb_control%recv_work_repl)) THEN
DO j = 1, SIZE(mixed_cdft%dlb_control%recv_work_repl)
IF (mixed_cdft%dlb_control%recv_work_repl(j)) THEN
DO i = 1, SIZE(mixed_cdft%dlb_control%recvbuff(j)%buffs)
index = (/LBOUND(mixed_cdft%dlb_control%recvbuff(j)%buffs(i)%cavity, 1), &
UBOUND(mixed_cdft%dlb_control%recvbuff(j)%buffs(i)%cavity, 1), &
LBOUND(mixed_cdft%dlb_control%recvbuff(j)%buffs(i)%cavity, 2), &
UBOUND(mixed_cdft%dlb_control%recvbuff(j)%buffs(i)%cavity, 2), &
LBOUND(mixed_cdft%dlb_control%recvbuff(j)%buffs(i)%cavity, 3), &
UBOUND(mixed_cdft%dlb_control%recvbuff(j)%buffs(i)%cavity, 3)/)
cdft_control_target%becke_control%cavity_mat(INDEX(1):INDEX(2), &
INDEX(3):INDEX(4), &
INDEX(5):INDEX(6)) = &
mixed_cdft%dlb_control%recvbuff(j)%buffs(i)%cavity
DEALLOCATE (mixed_cdft%dlb_control%recvbuff(j)%buffs(i)%cavity)
END DO
END IF
END DO
END IF
END IF
END IF
DO j = 1, SIZE(mixed_cdft%source_list)
DEALLOCATE (recvbuffer(j)%r3)
END DO
IF (calculate_forces) THEN
! Gradients of the weight function
DO j = 1, SIZE(mixed_cdft%source_list)
ALLOCATE (recvbuffer(j)%r4(3*natom, recvbuffer(j)%imap(1):recvbuffer(j)%imap(2), &
recvbuffer(j)%imap(3):recvbuffer(j)%imap(4), &
recvbuffer(j)%imap(5):recvbuffer(j)%imap(6)))
CALL force_env%para_env%irecv(msgout=recvbuffer(j)%r4, source=mixed_cdft%source_list(j), &
request=req_total(j))
END DO
DO i = 1, my_special_work
DO j = 1, SIZE(mixed_cdft%dest_list)
ind = j + (i - 1)*SIZE(mixed_cdft%dest_list) + SIZE(mixed_cdft%source_list)
IF (mixed_cdft%is_special) THEN
CALL force_env%para_env%isend(msgin=mixed_cdft%sendbuff(j)%gradients, &
dest=mixed_cdft%dest_list(j) + (i - 1)*force_env%para_env%num_pe/2, &
request=req_total(ind))
ELSE
CALL force_env%para_env%isend(msgin=cdft_control%group(1)%gradients, dest=mixed_cdft%dest_list(j), &
request=req_total(ind))
END IF
END DO
END DO
CALL mp_waitall(req_total)
IF (mixed_cdft%is_special) THEN
DO j = 1, SIZE(mixed_cdft%dest_list)
DEALLOCATE (mixed_cdft%sendbuff(j)%gradients)
END DO
DEALLOCATE (mixed_cdft%sendbuff)
ELSE
DEALLOCATE (cdft_control%group(1)%gradients)
END IF
ALLOCATE (cdft_control_target%group(1)%gradients(3*natom, bo(1, 1):bo(2, 1), &
bo(1, 2):bo(2, 2), lb_min:ub_max))
DO j = 1, SIZE(mixed_cdft%source_list)
cdft_control_target%group(1)%gradients(:, recvbuffer(j)%imap(1):recvbuffer(j)%imap(2), &
recvbuffer(j)%imap(3):recvbuffer(j)%imap(4), &
recvbuffer(j)%imap(5):recvbuffer(j)%imap(6)) = recvbuffer(j)%r4
DEALLOCATE (recvbuffer(j)%r4)
END DO
IF (mixed_cdft%dlb) THEN
IF (ANY(mixed_cdft%dlb_control%recv_work_repl)) THEN
DO j = 1, SIZE(mixed_cdft%dlb_control%recv_work_repl)
IF (mixed_cdft%dlb_control%recv_work_repl(j)) THEN
DO i = 1, SIZE(mixed_cdft%dlb_control%recvbuff(j)%buffs)
index = (/LBOUND(mixed_cdft%dlb_control%recvbuff(j)%buffs(i)%gradients, 2), &
UBOUND(mixed_cdft%dlb_control%recvbuff(j)%buffs(i)%gradients, 2), &
LBOUND(mixed_cdft%dlb_control%recvbuff(j)%buffs(i)%gradients, 3), &
UBOUND(mixed_cdft%dlb_control%recvbuff(j)%buffs(i)%gradients, 3), &
LBOUND(mixed_cdft%dlb_control%recvbuff(j)%buffs(i)%gradients, 4), &
UBOUND(mixed_cdft%dlb_control%recvbuff(j)%buffs(i)%gradients, 4)/)
cdft_control_target%group(1)%gradients(:, INDEX(1):INDEX(2), &
INDEX(3):INDEX(4), &
INDEX(5):INDEX(6)) = &
mixed_cdft%dlb_control%recvbuff(j)%buffs(i)%gradients
DEALLOCATE (mixed_cdft%dlb_control%recvbuff(j)%buffs(i)%gradients)
END DO
END IF
END DO
END IF
END IF
END IF
! Clean up remaining temporaries
IF (mixed_cdft%dlb) THEN
IF (ANY(mixed_cdft%dlb_control%recv_work_repl)) THEN
DO j = 1, SIZE(mixed_cdft%dlb_control%recv_work_repl)
IF (mixed_cdft%dlb_control%recv_work_repl(j)) THEN
IF (ASSOCIATED(mixed_cdft%dlb_control%recv_info(j)%target_list)) &
DEALLOCATE (mixed_cdft%dlb_control%recv_info(j)%target_list)
DEALLOCATE (mixed_cdft%dlb_control%recvbuff(j)%buffs)
END IF
END DO
DEALLOCATE (mixed_cdft%dlb_control%recv_info, mixed_cdft%dlb_control%recvbuff)
END IF
IF (ASSOCIATED(mixed_cdft%dlb_control%target_list)) &
DEALLOCATE (mixed_cdft%dlb_control%target_list)
DEALLOCATE (mixed_cdft%dlb_control%recv_work_repl)
END IF
DEALLOCATE (recvbuffer)
DEALLOCATE (req_total)
DEALLOCATE (lb)
DEALLOCATE (ub)
CALL timestop(handle2)
! Set some flags so the weight is not rebuilt during SCF
cdft_control_target%external_control = .TRUE.
cdft_control_target%need_pot = .FALSE.
cdft_control_target%transfer_pot = .FALSE.
! Store the bound indices for force calculation
IF (calculate_forces) THEN
cdft_control_target%becke_control%confine_bounds(2) = ub_max
cdft_control_target%becke_control%confine_bounds(1) = lb_min
END IF
CALL pw_scale(cdft_control_target%group(1)%weight, &
cdft_control_target%group(1)%weight%pw_grid%dvol)
! Set flags for ET coupling calculation
IF (mixed_env%do_mixed_et) THEN
IF (MODULO(mixed_cdft%sim_step, mixed_env%et_freq) == 0) THEN
dft_control%qs_control%cdft_control%do_et = .TRUE.
dft_control%qs_control%cdft_control%calculate_metric = mixed_cdft%calculate_metric
ELSE
dft_control%qs_control%cdft_control%do_et = .FALSE.
dft_control%qs_control%cdft_control%calculate_metric = .FALSE.
END IF
END IF
t2 = m_walltime()
IF (iounit > 0) THEN
WRITE (iounit, '(A)') ' '
WRITE (iounit, '(T2,A,F6.1,A)') 'MIXED_CDFT| Becke constraint built in ', t2 - t1, ' seconds'
WRITE (iounit, '(A)') ' '
END IF
CALL cp_print_key_finished_output(iounit, logger, force_env_section, &
"MIXED%MIXED_CDFT%PRINT%PROGRAM_RUN_INFO")
CALL timestop(handle)
END SUBROUTINE mixed_cdft_build_weight_parallel
! **************************************************************************************************
!> \brief Transfer CDFT weight/gradient between force_evals
!> \param force_env the force_env that holds the CDFT sub_force_envs
!> \param calculate_forces if forces should be computed
!> \param iforce_eval index of the currently active CDFT state
!> \par History
!> 01.2017 created [Nico Holmberg]
! **************************************************************************************************
SUBROUTINE mixed_cdft_transfer_weight(force_env, calculate_forces, iforce_eval)
TYPE(force_env_type), POINTER :: force_env
LOGICAL, INTENT(IN) :: calculate_forces
INTEGER, INTENT(IN) :: iforce_eval
CHARACTER(len=*), PARAMETER :: routineN = 'mixed_cdft_transfer_weight'
INTEGER :: bounds_of(8), handle, iatom, igroup, &
jforce_eval, nforce_eval
LOGICAL, SAVE :: first_call = .TRUE.
TYPE(cdft_control_type), POINTER :: cdft_control_source, cdft_control_target
TYPE(dft_control_type), POINTER :: dft_control_source, dft_control_target
TYPE(force_env_type), POINTER :: force_env_qs_source, force_env_qs_target
TYPE(mixed_cdft_type), POINTER :: mixed_cdft
TYPE(mixed_environment_type), POINTER :: mixed_env
TYPE(pw_env_type), POINTER :: pw_env_source, pw_env_target
TYPE(pw_pool_type), POINTER :: auxbas_pw_pool_source, &
auxbas_pw_pool_target
TYPE(qs_environment_type), POINTER :: qs_env_source, qs_env_target
NULLIFY (mixed_cdft, dft_control_source, dft_control_target, force_env_qs_source, &
force_env_qs_target, pw_env_source, pw_env_target, auxbas_pw_pool_source, &
auxbas_pw_pool_target, qs_env_source, qs_env_target, mixed_env, &
cdft_control_source, cdft_control_target)
mixed_env => force_env%mixed_env
CALL get_mixed_env(mixed_env, cdft_control=mixed_cdft)
CALL timeset(routineN, handle)
IF (iforce_eval == 1) THEN
jforce_eval = 1
ELSE
jforce_eval = iforce_eval - 1
END IF
nforce_eval = SIZE(force_env%sub_force_env)
SELECT CASE (force_env%sub_force_env(jforce_eval)%force_env%in_use)
CASE (use_qs_force, use_qmmm)
force_env_qs_source => force_env%sub_force_env(jforce_eval)%force_env
force_env_qs_target => force_env%sub_force_env(iforce_eval)%force_env
CASE DEFAULT
CPASSERT(.FALSE.)
END SELECT
IF (.NOT. force_env%mixed_env%do_mixed_qmmm_cdft) THEN
CALL force_env_get(force_env=force_env_qs_source, &
qs_env=qs_env_source)
CALL force_env_get(force_env=force_env_qs_target, &
qs_env=qs_env_target)
ELSE
qs_env_source => force_env_qs_source%qmmm_env%qs_env
qs_env_target => force_env_qs_target%qmmm_env%qs_env
END IF
IF (iforce_eval == 1) THEN
! The first force_eval builds the weight function and gradients in qs_cdft_methods.F
! Set some flags so the constraint is saved if the constraint definitions are identical in all CDFT states
CALL get_qs_env(qs_env_source, dft_control=dft_control_source)
cdft_control_source => dft_control_source%qs_control%cdft_control
cdft_control_source%external_control = .FALSE.
cdft_control_source%need_pot = .TRUE.
IF (mixed_cdft%identical_constraints) THEN
cdft_control_source%transfer_pot = .TRUE.
ELSE
cdft_control_source%transfer_pot = .FALSE.
END IF
mixed_cdft%sim_step = mixed_cdft%sim_step + 1
ELSE
! Transfer the constraint from the ith force_eval to the i+1th
CALL get_qs_env(qs_env_source, dft_control=dft_control_source, &
pw_env=pw_env_source)
CALL pw_env_get(pw_env_source, auxbas_pw_pool=auxbas_pw_pool_source)
cdft_control_source => dft_control_source%qs_control%cdft_control
CALL get_qs_env(qs_env_target, dft_control=dft_control_target, &
pw_env=pw_env_target)
CALL pw_env_get(pw_env_target, auxbas_pw_pool=auxbas_pw_pool_target)
cdft_control_target => dft_control_target%qs_control%cdft_control
! The constraint can be transferred only when the constraint defitions are identical in all CDFT states
IF (mixed_cdft%identical_constraints) THEN
! Weight function
DO igroup = 1, SIZE(cdft_control_target%group)
ALLOCATE (cdft_control_target%group(igroup)%weight)
CALL auxbas_pw_pool_target%create_pw(cdft_control_target%group(igroup)%weight)
! We have ensured that the grids are consistent => no danger in using explicit copy
CALL pw_copy(cdft_control_source%group(igroup)%weight, cdft_control_target%group(igroup)%weight)
CALL auxbas_pw_pool_source%give_back_pw(cdft_control_source%group(igroup)%weight)
DEALLOCATE (cdft_control_source%group(igroup)%weight)
END DO
! Cavity
IF (cdft_control_source%type == outer_scf_becke_constraint) THEN
IF (cdft_control_source%becke_control%cavity_confine) THEN
CALL auxbas_pw_pool_target%create_pw(cdft_control_target%becke_control%cavity)
CALL pw_copy(cdft_control_source%becke_control%cavity, cdft_control_target%becke_control%cavity)
CALL auxbas_pw_pool_source%give_back_pw(cdft_control_source%becke_control%cavity)
END IF
END IF
! Gradients
IF (calculate_forces) THEN
DO igroup = 1, SIZE(cdft_control_source%group)
bounds_of = (/LBOUND(cdft_control_source%group(igroup)%gradients, 1), &
UBOUND(cdft_control_source%group(igroup)%gradients, 1), &
LBOUND(cdft_control_source%group(igroup)%gradients, 2), &
UBOUND(cdft_control_source%group(igroup)%gradients, 2), &
LBOUND(cdft_control_source%group(igroup)%gradients, 3), &
UBOUND(cdft_control_source%group(igroup)%gradients, 3), &
LBOUND(cdft_control_source%group(igroup)%gradients, 4), &
UBOUND(cdft_control_source%group(igroup)%gradients, 4)/)
ALLOCATE (cdft_control_target%group(igroup)% &
gradients(bounds_of(1):bounds_of(2), bounds_of(3):bounds_of(4), &
bounds_of(5):bounds_of(6), bounds_of(7):bounds_of(8)))
cdft_control_target%group(igroup)%gradients = cdft_control_source%group(igroup)%gradients
DEALLOCATE (cdft_control_source%group(igroup)%gradients)
END DO
END IF
! Atomic weight functions needed for CDFT charges
IF (cdft_control_source%atomic_charges) THEN
IF (.NOT. ASSOCIATED(cdft_control_target%charge)) &
ALLOCATE (cdft_control_target%charge(cdft_control_target%natoms))
DO iatom = 1, cdft_control_target%natoms
CALL auxbas_pw_pool_target%create_pw(cdft_control_target%charge(iatom))
CALL pw_copy(cdft_control_source%charge(iatom), cdft_control_target%charge(iatom))
CALL auxbas_pw_pool_source%give_back_pw(cdft_control_source%charge(iatom))
END DO
END IF
! Set some flags so the weight is not rebuilt during SCF
cdft_control_target%external_control = .FALSE.
cdft_control_target%need_pot = .FALSE.
! For states i+1 < nforce_eval, prevent deallocation of constraint
IF (iforce_eval == nforce_eval) THEN
cdft_control_target%transfer_pot = .FALSE.
ELSE
cdft_control_target%transfer_pot = .TRUE.
END IF
cdft_control_target%first_iteration = .FALSE.
ELSE
! Force rebuild of constraint and dont save it
cdft_control_target%external_control = .FALSE.
cdft_control_target%need_pot = .TRUE.