forked from alexei-matveev/paragauss-gpl
-
Notifications
You must be signed in to change notification settings - Fork 0
/
main_scf.f90
1634 lines (1430 loc) · 56.3 KB
/
main_scf.f90
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
!
! ParaGauss, a program package for high-performance computations of
! molecular systems
!
! Copyright (C) 2014 T. Belling, T. Grauschopf, S. Krüger,
! F. Nörtemann, M. Staufer, M. Mayer, V. A. Nasluzov, U. Birkenheuer,
! A. Hu, A. V. Matveev, A. V. Shor, M. S. K. Fuchs-Rohr, K. M. Neyman,
! D. I. Ganyushin, T. Kerdcharoen, A. Woiterski, A. B. Gordienko,
! S. Majumder, M. H. i Rotllant, R. Ramakrishnan, G. Dixit,
! A. Nikodem, T. Soini, M. Roderus, N. Rösch
!
! This program is free software; you can redistribute it and/or modify
! it under the terms of the GNU General Public License version 2 as
! published by the Free Software Foundation [1].
!
! This program is distributed in the hope that it will be useful, but
! WITHOUT ANY WARRANTY; without even the implied warranty of
! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
! General Public License for more details.
!
! [1] http://www.gnu.org/licenses/gpl-2.0.html
!
! Please see the accompanying LICENSE file for further information.
!
!=====================================================================
! Public interface of module
!=====================================================================
subroutine main_scf()
!
! Purpose: Main Level of the SCF part. Executed on all workers.
!
! Subroutine called by: main_master
!
! Author: Folke Noertemann
! Date: 10/95
!
!===================================================================
! End of Public interface of module
!===================================================================
!-------------------------------------------------------------------
! Modifications
!-------------------------------------------------------------------
!
! Author: UB
! Date: 27/5/97
! Description: Introduction of the "USE_MODEL_DENSITY" option
!
! Modification (Please copy before editing)
! Author: ...
! Date: ...
! Description: ...
!
! Modification (Please copy before editing)
! Author: AS
! Date: 7/98
! Description: pvm -> comm
!
! Modification (Please copy before editing)
! Author: ...
! Date: ...
! Description: ...
!
!-------------------------------------------------------------------
!------------ Modules used -----------------------------------------
# include "def.h"
use type_module, only: i4_kind, r8_kind ! type specification parameters
use interfaces, only: chargefit
use paragauss, only: toggle_legacy_mode
use prescf_module, only: prescf_init, prescf_finalize
use time_module, only: start_timer, stop_timer ! timing routines
use timer_module, only: timer_grid_setup, timer_scf, &
timer_scf_cycle, timer_scf_ham, timer_scf_preparations, &
timer_scf_reoc, timer_scf_xc, timer_grid_small_to_large, &
timer_print_scfcycle, timer_print_scf
use occupation_module, only: fixed_hole, n_occo, reoccup,&
print_occ_num, occupation_2d_correct, &
occupation_print_spectrum, occupation_get_holes
use fermi_module, only: fermi_reoccup, fermi_level_broad, fermi_get_entropy, &
fermi_write_input, fermi_read_scfcontrol
use occupied_levels_module, only: sndrcv_eigvec_occ
use options_module, only: xcmode_model_density, xcmode_extended_mda, &
xcmode_numeric_exch, recover_fitcoeff, recover_scfstate, recover_ksmatrix, &
xcmode_exchange_fit, recover_eigenvec, recover_fragment, options_save_ksmatrix, &
options_save_scfstate, options_save_interval, options_save_eigenvec_all, &
options_save_eigenvec, options_save_densmatrix, options_save_fitcoeff, &
options_xcmode, options_recover, lvshift_mixing, options_spin_orbit, &
recover_nothing, options_save_as_fragment, options_reset_scfcycle
use output_module, only: output_data_read, output_data_saved, output_timing_scfloops, &
output_scfloops, output_timing_scf, output_timing_detailedscf, &
output_spectrum, output_reoccup, output_each_eigendata, output_main_scf, &
output_eigendata, output_densmat, output_chargefit
use iounitadmin_module, only: get_iounit, return_iounit, output_unit, stdout_unit, &
write_to_trace_unit, write_to_output_units, openget_iounit, returnclose_iounit
use grid_module, only: grid_main, grid_close
use xc_cntrl, only: xc_is_on=>is_on, xc_ANY
use xc_hamiltonian, only: xc_setup, xc_build, &
xc_hamiltonian_store, xc_hamiltonian_recover, xc_close, s_average
use xcfit_hamiltonian, only: xcfit_setup, xcfit_build, xcfit_close
use xcmda_hamiltonian, only: xcmda_setup, xcmda_coeff_store, &
xcmda_build, xcmda_coeff_recover, xcmda_close, &
rho_exact, comp_exact
use convergence_module, only: convergence_abort_calculation, convergence, &
convergence_check_density_dev, convergence_max_iterations, &
convergence_state_store, convergence_check_coulomb_dev, &
convergence_check_coeff_dev, convergence_state_recover, &
convergence_clear_buffers, convergence_setup, convergence_resize_buffers, &
convergence_put_energy, convergence_put_density_dev, convergence_put_coeff_dev, &
convergence_put_coulomb_dev, convergence_shutdown
use eigen_data_module, only: print_eigendata, eigvec_write, eigen_data_dump
use fit_coeff_module, only: print_coeff_charge, &
fit_coeff_store, fit_coeff_recover, fit_coeff_sndrcv, &
fit_coeff_normalize, ICHFIT
#ifdef WITH_CORE_DENS
use fit_coeff_module, only: write_coeff_core
#endif
use density_data_module, only: print_densmat, gendensmat_occ, &
save_densmat, densmat
use overlap_module, only: overlap
use filename_module, only: outfile, recfile
use hamiltonian_module, only: ham_tot, hamiltonian_store, hamiltonian_recover
use eigen_data_module, only: eigen_data_solve, eigen_hole_shutdown, &
build_lvsft_hamiltonian, make_level_shift
use diis_fock_module, only: diis_fock_matrix, fixed_fock_matrix, &
diis_fock_module_close, diis_on
use energy_calc_module, only: write_energies, get_energy
use operations_module, only: operations_geo_opt, &
operations_core_density, operations_solvation_effect, operations_qm_mm
use symmetry_data_module, only: symmetry_data_n_spin
use ham_calc_module, only: ham_calc_main
use mixing_module, only: mixing_state_store, mixing_state_recover, &
mixing_close, mixing_write_input, mixing_reset_buffers, &
mixing_beta_ch, mixing_beta_sp, mixing_beta_xc, level_shift, &
start_after_cycle
#ifdef WITH_SCFCONTROL
use mixing_module, only: mixing_read_scfcontrol
use convergence_module, only: convergence_write, convergence_read_scfcontrol
use options_module, only : options_write_input, options_read_scfcontrol
#endif
use population_module, only: population_mulliken, population_spor_mulliken
use comm, only: comm_parallel, comm_same, comm_bcast, comm_rank
use pairs_module, only: n_pairs
use readwriteblocked_module, only: readwriteblocked_tapehandle, &
readwriteblocked_stopread, readwriteblocked_startread, &
readwriteblocked_startwrite, readwriteblocked_read, readwriteblocked_stopwrite, &
readwriteblocked_write
use solv_electrostat_module, only: build_solv_ham, &
dealloc_ham_solv, solv_energy_el, calc_q_e
use solv_cavity_module, only: sol_start_cycle, n_Q_update
use solv_charge_mixing_module, only: mix_charges, Qs_mix, Q_dealloc
use induced_dipoles_module, only: calc_Pol_centers, free_field_arrays, n_update
use calc_id_module, only: build_Pol_ham, calc_id_energy, dealloc_Pol_ham
#ifdef WITH_EFP
use induced_dipoles_module, only: do_pol_pcm
use efp_module, only: n_efp, qm_fixed, print_id
use efp_polar_module, only: deallocate_Efield
use efp_solv_module, only: allocate_V_and_Q_id, calc_V_and_Q_id, allocate_E_cav, do_solv
use efp_solv_module, only: deallocate_E_cav, deallocate_V_id
#endif
use s2_expect, only: s2_calc
implicit none
! *** end of interface ***
!------------ Declaration of local variables --------------------
integer (i4_kind) :: xcmode, recover_mode, loop
integer (i4_kind) :: first_loop, trace_width(14)
real (r8_kind) :: tot_en, spin_diff, time_for_last_cycle, entropy
real (r8_kind) :: coeff_dev, coulomb_dev, density_dev
character (len=16) :: trace_format(15), trace_line*101
integer, parameter :: NCHAR = 18 ! string length
character (len=NCHAR), parameter :: fit_file = 'saved_fitcoeff.dat'
character (len=NCHAR), parameter :: scf_file = 'saved_scfstate.dat'
character (len=NCHAR), parameter :: ham_file = 'saved_ksmatrix.dat'
character (len=NCHAR), parameter :: eig_file = 'saved_eigenvec.dat'
character (len=NCHAR), parameter :: frg_file = 'saved_fragment.dat'
! One of the above, including .dat extension:
character (len=NCHAR) :: rec_file
logical :: scf_conv ! signifies SCF convergence, broadcasted to slaves
logical :: check_density_dev
logical :: store_now, recover, etot_recovered
logical :: new_trace_header
logical :: update_scfcontrol = .FALSE.
integer (i4_kind), save :: scf_count = 0 ! counts calls to main_scf, i.e. geo_loops
!ifdef WITH_EFP
logical :: do_update
!endif
!
! As of now main_scf() runs on all workers. If you are thinking
! about adding code to be executed on all workers, consider adding
! code/calls outside of the legacy mode blocks or to other subs
! executed in parallel:
!
! - prescf_init()
! - grid_main()
! - ham_calc_main()
! - prescf_finalize()
! - ...
!
call say ("start")
! We are called once again, remember that:
scf_count = scf_count + 1
! Write to trace unit to show progress of calculation
call write_to_trace_unit ("Entering SCF part")
call start_timer (timer_scf)
call start_timer (timer_scf_preparations)
! Check status of recover files (before entering prescf and
! xc_setup_main)
call init_recover (recover, recover_mode, rec_file) ! all intent (out)
! Preparatory tasks for scf-cycle
call say ("prescf_init")
!
! This sub is executed on all workers, put the code that is common
! to all workers inside (or rather call from there):
!
call prescf_init()
xcmode = options_xcmode()
! Prepare calculation of exchange-correlation hamiltonian:
if (xc_is_on (xc_ANY)) then
! Set up grid
call say ("grid setup")
call start_timer (timer_grid_setup)
call say ("call grid_main()")
call grid_main (.false.)
call say ("done grid_main()")
call stop_timer (timer_grid_setup)
! Call xc setup:
call say ("call xc_setup()")
select case (xcmode)
case (xcmode_model_density, xcmode_extended_mda)
call xcmda_setup ()
case (xcmode_numeric_exch)
call xc_setup ()
case (xcmode_exchange_fit)
call xcfit_setup ()
end select
call say ("done xc_setup()")
endif
if (output_chargefit) then
call say ("coeff_charge is written to file")
call print_coeff_charge (0)
endif
call say ("other preparations")
! Reset convergence test variables and allocated the convergence
! buffers
call convergence_setup()
#ifdef WITH_SCFCONTROL
! Prepare file with steering parameters read at each cycle
if (scf_count .eq. 1) then
call write_scfcontrol()
else
call read_scfcontrol()
endif
#endif /* ifdef WITH_SCFCONTROL */
! It was previousely done every SCF iteration in
! convergence_read_scfcontrol() called from read/write_scfcontrol
call convergence_resize_buffers()
call stop_timer (timer_scf_preparations)
call say ("preparations done")
call write_to_trace_unit ("SCF preparations done")
! Here the SCF-Cycle starts
call write_trace_header()
time_for_last_cycle = 0.0_r8_kind
! FIXME: what about tot_en?
loop = 0
etot_recovered = .false.
ASSERT(comm_same(recover))
if (recover) then
!
! Loads loop, etot_recovered and tot_en. When calling
! hamiltonian_recover() the Fock matrix in ham_tot is allocated
! and filled with data. Master seems to send the orders and data
! to slaves from there, the correspondig receives are initiated
! from main_slave():
!
call do_recover (recover_mode, n_pairs, loop, tot_en, eof=etot_recovered)
! loop is incremented at the top of the scf_cycle:
loop = loop - 1
endif
! This is compared against base-1 loop:
first_loop = loop + 1
ASSERT(comm_same(fixed_hole))
if (fixed_hole) then
call say ("set up hole information in eigen_data_module")
call occupation_get_holes()
endif
scf_cycle: do
loop = loop + 1
call start_timer (timer_scf_cycle)
! Saving state every iteration may require much IO, do it maybe
! every few iterations:
store_now = store_now_p (loop, options_save_interval())
if (output_timing_scfloops .or. output_scfloops) then
call write_to_output_units (" ")
call write_to_output_units ("####### Loop: ", loop)
call write_to_output_units (" ")
endif
! Check recover options and fork accordingly
ASSERT(comm_same(loop))
ASSERT(comm_same(first_loop))
if (loop == first_loop .and. recover) then
ASSERT(comm_same(recover_mode))
select case (recover_mode)
case (recover_fitcoeff); goto 1000
case (recover_scfstate); goto 1000
case (recover_ksmatrix); goto 2000
case (recover_eigenvec); goto 3000
case (recover_fragment); goto 3000
end select
endif
! Calculation of exchange-correlation hamiltonian:
if (xc_is_on (xc_ANY)) then
ASSERT(comm_same(xcmode))
select case (xcmode)
case (xcmode_model_density, xcmode_extended_mda)
ABORT("verify SPMD")
call say ("fit_coeff_sndrcv")
call fit_coeff_sndrcv (ICHFIT)
if ((rho_exact .or. comp_exact) .and. loop /= first_loop) then
call sndrcv_eigvec_occ()
endif
call start_timer (timer_scf_xc)
call say ("xcmda_build_ham")
call xcmda_build ((loop == first_loop))
call stop_timer (timer_scf_xc)
case (xcmode_numeric_exch)
if (loop /= first_loop) then
call start_timer (timer_scf_xc)
call say ("xc_build")
! Call with relative cycle number:
call xc_build (loop - first_loop + 1)
call stop_timer (timer_scf_xc)
endif
case (xcmode_exchange_fit)
if (loop /= first_loop) then
call start_timer (timer_scf_xc)
call say ("xcfit_build")
call xcfit_build (loop)
call stop_timer (timer_scf_xc)
endif
end select
endif
! The code inside this DO block is executed on master only. The
! slaves are spinning in main_slave() and waiting for orders
! during that time:
do while (toggle_legacy_mode ())
! Master-only context here ...
if (calc_Pol_centers() .and. loop >= first_loop+1) then
if (loop == first_loop+1) then
do_update = .true.
else
do_update = (mod (loop - (first_loop + 1), n_update) == 0)
end if
#ifdef WITH_EFP
if (operations_solvation_effect) then
if (do_pol_pcm .and. n_efp > 0) then
call allocate_V_and_Q_id()
call allocate_E_cav()
end if
! This is a variable in a efp_solv_module:
do_solv = (loop > first_loop + sol_start_cycle)
end if
call build_Pol_ham (print_id, do_update)
#else
call build_Pol_ham (do_update)
#endif
end if
enddo ! while (toggle_legacy_mode())
! Parallel context from here on ...
! Calculating solvation hamiltonian
if (operations_solvation_effect .and. &
loop >= first_loop+sol_start_cycle) then
if (mod (loop - (first_loop + sol_start_cycle), n_Q_update) == 0) then
call calc_Q_e()
end if
#ifdef WITH_EFP
if (do_pol_pcm .and. n_efp > 0) then
call calc_V_and_Q_id()
end if
#endif
call build_solv_ham (loop, first_loop + sol_start_cycle)
endif
if (output_densmat) then
call say ("densmat is written to file")
call print_densmat (loop)
endif
! Entry point for "read_fitcoeff" and "read_scfstate":
1000 continue
! Save fit coefficients and current SCF state if required
if (options_save_fitcoeff()) then
call do_fitcoeff_store (store_now) ! reads loop
endif
! Save current SCF state if required
if (options_save_scfstate()) then
call do_scfstate_store (store_now) ! reads loop
endif
! Calculation of other (not exchange correlation) parts of
! hamiltonian and summation to total hamiltonian and calculation
! of energies
call say ("ham_calc_main")
call start_timer (timer_scf_ham)
!
! Build Coulomb and Fock matrices. This also call reset_ham()
! which allocates and initializes Fock matrix in ham_tot. Some
! recover options jump over this call.
!
call ham_calc_main (loop - first_loop + 1) ! relative cycle number
call stop_timer (timer_scf_ham)
if (operations_solvation_effect .and. &
loop >= first_loop + sol_start_cycle) then
! FIXME: slaves did it earlier, see build_solv_ham(), but
! since it is idempotent we just call it again:
call dealloc_ham_solv () ! no comm, idempotent
! This one does not do much on slaves:
call solv_energy_el () ! no comm
endif
if (calc_Pol_centers() .and. loop >= first_loop+1) then
do while (toggle_legacy_mode ())
! Master-only context here ...
call dealloc_Pol_ham()
call calc_id_energy()
enddo
end if
! Write calculated energies. Slaves could always print that to
! stdout, however the values are partial there.
call say ("write_energies")
if (output_scfloops .or. output_main_scf) then
call write_energies (output_unit)
if (comm_rank() == 0) then
call write_energies (stdout_unit)
endif
endif
if (etot_recovered) then
if (output_scfloops .or. output_main_scf) then
500 format(2X, 'e_sum = ', F25.12, ' (recovered)')
write (output_unit, 500) tot_en
if (comm_rank() == 0) then
write (stdout_unit, 500) tot_en
endif
endif
else
call get_energy (tot=tot_en)
endif
! And pass energy to the convergence check module
call convergence_put_energy (tot_en)
if (store_now) then
! Write tot_en to a file:
call store_total_energy()
endif
! Entry point for "read_ksmatrix":
2000 continue
! Save Kohn-Sham matrix if required. NOOP on slaves:
if (options_save_ksmatrix()) then
call do_ksmatrix_store (store_now) ! reads loop
endif
#ifdef WITH_SCFCONTROL
#warning "WITH_SCFCONTROL=1 is not useful for batch runs"
! Read file with steering parameters first and modify the length
! of convergence buffers (if required) keeping a copy of the
! original data. Caution: reloads options_save_interval(),
! options_save_fitcoeff(). If you want every worker to know
! updated parameters, you either have to let everyone read them,
! or broadcast them.
call read_scfcontrol()
#endif /* ifdef WITH_SCFCONTROL */
! Check convergence and exit loop in case of convergence
call say ("convergence is checked")
! Check SCF convergence. FIXME: maybe make sure all workers have
! the data to determine convergence on their own?
scf_conv = (convergence() .or. convergence_max_iterations (loop)) &
.and. loop /= first_loop
call comm_bcast (scf_conv)
if (scf_conv) then
call stop_timer (timer_scf_cycle)
call timer_grid_small_to_large()
! Does nothing when the slaves have no output:
if (output_timing_scfloops) then
call timer_print_scfcycle()
endif
! Uses global variables loop, tot_en and old values of the
! save_data options. NOOP on slaves.
call do_final_store (n_vir=n_pairs)
call say ("exit scf_cycle")
! Let everyone exit the SCF block. This is the only exit
! point:
exit scf_cycle
endif
! Now update the save_data options, deallocate the temporarily
! kept convergence buffers (as far as necessary), reset the
! mixing buffers (if required) and update the scf_control file
! accordingly
store_now = store_now_p (loop, options_save_interval())
call convergence_clear_buffers()
call mixing_reset_buffers (loop, update_scfcontrol)
! Note that update_scfcontrol is also set by read_scfcontrol
#ifdef WITH_SCFCONTROL
if (update_scfcontrol) call write_scfcontrol()
#endif /* ifdef WITH_SCFCONTROL */
!
! In diis_fock_matrix() the total hamiltonian could be updated in
! a DIIS routine to get a convergence acceleration. Only ham_tot
! will be changed by the routine. It could be switched on by
! diis_on = true in the DIIS input namelist. Parameters need to
! be adapted; their best choice depend on other mixing strategies
! like Perturbation Theory and Fermi broadening.
!
if (diis_on) then
ASSERT (.not.options_spin_orbit)
ASSERT (allocated(overlap))
ASSERT (allocated(ham_tot))
ASSERT (allocated(densmat))
call diis_fock_matrix (ham_tot, overlap, densmat, loop)
call fixed_fock_matrix (ham_tot, loop)
endif
! Eigensolver to determine orbital coefficients
call say ("eigen_data_solve")
!
! Perform Hamiltonian diagonalization, see eigen_data_module.
! FIXME: pass arguments explicitly.
!
call eigen_data_solve()
! Write eigendata to temporary file "eigenmist" for debugging
if (output_each_eigendata) then
call say ("print_eigendata")
call print_eigendata()
endif
! Reoccupation of orbitals
call start_timer (timer_scf_reoc)
if (fermi_level_broad) then
call say ("fermi_reoccup")
! FIXME: The input differes on master and slaves so that
! Newton iteration fails:
if (comm_rank() == 0) then
call fermi_reoccup()
endif
else
call say ("reoccup")
call reoccup()
endif
call say ("occupation_2d_correct")
if (.not. options_spin_orbit) then
! Does no communication. Rewrites occ_num(), n_occo() in
! place:
call occupation_2d_correct()
endif
call stop_timer (timer_scf_reoc)
! Call on master after eigen_data_solve() and reocupp(). FIXME:
! The procedure does no communication, but fiddles with global
! variables that may not be available on slaves. Verify this next
! time the ABORT() fires.
if (lvshift_mixing) then
ABORT("verify SPMD")
call build_lvsft_hamiltonian &
(n_occo, level_shift, loop .gt. START_AFTER_CYCLE)
endif
if (output_spectrum) then
call say ("print out spectrum")
call occupation_print_spectrum()
endif
if (output_reoccup .or. output_each_eigendata) then
call say ("print_occ_num")
call print_occ_num (loop)
endif
! Entry point for "read_eigenvec":
3000 continue
! Save eigenstates and occupation if required. NOOP on slaves.
if (options_save_eigenvec()) then
call do_eigenvec_store (store_now, n_vir=n_pairs) ! reads loop
endif
! Setup of density matrix (on each slave) and fit charge
! fitfunction coefficients
call say ("send_eigvec_occ")
call sndrcv_eigvec_occ()
call say ("gendensmat_occ")
! The function here seems to just return the input parameter and
! thus should give the same value on all workers, modulo bugs, as
! always:
check_density_dev = convergence_check_density_dev() &
.and. .not. (recover .and. loop == first_loop)
ASSERT(comm_same(check_density_dev))
if (check_density_dev) then
call gendensmat_occ (density_dev)
else
call gendensmat_occ()
if (recover .and. loop == first_loop) then
density_dev = 0.0_r8_kind
else
density_dev = huge (0.0_r8_kind)
endif
endif
call convergence_put_density_dev (density_dev)
!
! Moved from gendensmat_occ() here:
!
if (output_densmat) then
call print_densmat (loop)
endif
call say ("chargefit")
!
! One of the first thigs it does is sending a message to the
! slaves telling them to also call chargefit() with dummy
! arguments that are assumed to be unused on slaves:
!
call chargefit (loop, coeff_dev, coulomb_dev)
call convergence_put_coeff_dev (coeff_dev)
call convergence_put_coulomb_dev (coulomb_dev)
if (output_chargefit) then
call say ("coeff_charge is written to file")
call print_coeff_charge (loop)
endif
! Calculate and print timing summary for SCF cycle
call say ("timing")
call stop_timer (timer_scf_cycle)
time_for_last_cycle = timer_scf_cycle%real_timediff
call timer_grid_small_to_large()
if (output_timing_scfloops) call timer_print_scfcycle()
! Write to trace unit to show progress of calculation
if (new_trace_header) then
call write_trace_header()
endif
call write_trace_item()
!
! In the next loop the energy is computed, not recovered:
!
etot_recovered = .false.
call say ("cycle done")
enddo scf_cycle
! This checks convergence of energy AND charge coefficients AND
! density matrix AND Coulomb (?). FIXME: This function seems to
! complete without crash but returns junk on slave workers:
scf_conv = convergence()
call comm_bcast (scf_conv)
if (scf_conv) then
call write_to_output_units ("")
call write_to_output_units ("####### convergence was reached in cycle ", inte=loop)
call write_to_output_units ("")
call write_to_trace_unit ("MAIN_SCF: convergence was reached in cycle ", inte=loop)
else
call write_to_output_units ("")
call write_to_output_units ("####### aborting at maximal number of cycles: ", inte=loop)
call write_to_output_units ("")
call write_to_trace_unit ("MAIN_SCF: aborting at maximal number of cycles: ", inte=loop)
endif
! Compute expectation values of S2. Does some output. FIXME: looks
! like we have a serial O(N^3) step here:
ASSERT(allocated(n_occo))
if (size (n_occo, 1) == 2) then ! unrestricted
call s2_calc() ! from s2_expect.f90
endif
#ifdef WITH_CORE_DENS
if (operations_core_density) then
ABORT("verify SPMD")
call write_coeff_core (loop)
end if
#endif
!
! The diis_fock_matrix() routine will not be used anymore in this
! SCF, so we deallocate the storage matrices. The
! diis_charge_coeff() routine is treated the same way. Make sure
! slaves do no harm when executing it! FIXME: maybe move it into
! prescf_finalize()?
!
call diis_fock_module_close()
! Clean up the XC part:
if (xc_is_on (xc_ANY)) then
call say ("call grid_close()")
call grid_close (.true.)
call say ("done grid_close()")
call say ("call xc_close()")
select case (xcmode)
case (xcmode_model_density, xcmode_extended_mda)
call xcmda_close ()
case (xcmode_numeric_exch)
call xc_close ()
case (xcmode_exchange_fit)
call xcfit_close ()
end select
call say ("done xc_close()")
endif
if (operations_solvation_effect) then
call calc_Q_e ()
call solv_energy_el () ! no comm
call Q_dealloc () ! no comm, idempotent
endif
if (calc_Pol_centers()) then
call calc_id_energy()
#ifdef WITH_EFP
if (.not. qm_fixed) call free_field_arrays()
#else
call free_field_arrays()
#endif
#ifdef WITH_EFP
if (n_efp > 1 .and. .not. qm_fixed) then
call deallocate_Efield()
end if
if (n_efp > 0 .and. operations_solvation_effect .and. do_pol_pcm) then
call deallocate_V_id()
call deallocate_E_cav()
end if
#endif
end if
! Write final calculated energies
call say ("write final calculated energies")
call write_energies (output_unit)
if (comm_rank() == 0) then
if (output_main_scf) call write_energies (stdout_unit)
endif
call get_energy (tot=tot_en)
! Write to trace unit to show progress of calculation
call write_to_trace_unit ("MAIN_SCF: final total energy: ", real=tot_en)
call say ("population analysis in progress ...")
if (options_spin_orbit) then
call population_spor_mulliken()
else
call population_mulliken()
endif
! Write eigendata to temporary file "eigvec.dat" for debugging
if (output_eigendata) then
call say ("print_eigendata")
call print_eigendata()
!
! Dumps all of hamiltonian, overlap, eigenvectors and eigenvalues
! into the named file in the output directory:
!
call eigen_data_dump ("eigendata.txt")
if (options_spin_orbit) then
DPRINT 'main_scf: call eigvec_write()'
call eigvec_write()
DPRINT 'main_scf: .'
endif
endif
if (lvshift_mixing) then
make_level_shift=.false.
endif
if (options_save_densmatrix()) call save_densmat !!!!!!!!!!!!AS
call say ("call prescf_finalize()")
call prescf_finalize()
call say ("done prescf_finalize()")
call say ("call mixing_close()")
call mixing_close (loop, xcmode == xcmode_exchange_fit)
call say ("done mixing_close()")
if (fixed_hole) then
call say ("shutdown hole information in eigen_data_module")
call eigen_hole_shutdown()
endif
! Print occupation of orbitals:
call say ("call occupation_print_spectrum()")
call occupation_print_spectrum ()
call say ("done occupation_print_spectrum()")
call stop_timer (timer_scf)
! Print timing summary for scf part:
if (output_timing_scf .or. output_timing_detailedscf) then
call say ("timer_print_scf()")
call timer_print_scf()
endif
if (convergence_abort_calculation()) &
call error_handler (&
"MAIN_SCF: aborting calculation because no convergence was reached")
! Deallocates the convergence test buffers:
call say ("convergence_shutdown")
call convergence_shutdown()
call say ("end")
contains
function store_now_p (loop, save_interval) result (yes)
implicit none
integer, intent (in) :: loop, save_interval
logical :: yes
! *** end of interface ***
if (save_interval /= 0) then
yes = mod (loop, save_interval) == 0
else
yes = .false.
endif
end function store_now_p
#ifdef WITH_SCFCONTROL
subroutine write_scfcontrol()
!
! Purpose: write variables that can be changed at each cycle to
! file. This function when executed in parallel, does not
! necessarily try writing to the same file. The outfile() returns
! a path to a scratch directory for slaves, different for each
! worker.
!
! To avoid warning messages due to spurious modifications of the
! scf control data caused by the finite accuracy achievable with
! formatted write and re-read, the scf control parametesr are
! re-read immediatelly after writing (UB, 8/97)
!
use iounitadmin_module, only: openget_iounit, returnclose_iounit
implicit none
!** End of interface ***************************************
integer :: unit
! not really an output file, but happens to reside in output directory:
unit = openget_iounit (trim (outfile ("scfcontrol")), status='replace')
call convergence_write (unit, full=.true.)
call options_write_input (unit, scfcontrol=.true.)
call mixing_write_input (unit, scfcontrol=.true.)
call fermi_write_input (unit, scfcontrol=.true.)
call returnclose_iounit (unit, status='keep')
call read_scfcontrol (warning=.false.)
update_scfcontrol = .false.
end subroutine write_scfcontrol
subroutine read_scfcontrol (warning)
!
! Purpose: read variables that can be changed at each cycle from
! file. FIXME: as is it looks like each worker would be reading
! different files as outfile() points to scratch for all but the
! rank-0 worker.
!
use input_module, only: input_open, input_close
logical, optional, intent (in) :: warning
!** End of interface ***************************************
! not really an output file, but happens to reside in output directory:
call input_open (trim (outfile ("scfcontrol")))
call convergence_read_scfcontrol (warning, new_trace_header)
call options_read_scfcontrol (warning)
call mixing_read_scfcontrol (warning)
call fermi_read_scfcontrol (warning, new_trace_header, update_scfcontrol)
call input_close()
end subroutine read_scfcontrol
#endif /* ifdef WITH_SCFCONTROL */
subroutine init_recover (recover, recover_mode, rec_file)
!
! Check recover files and set recover options accordingly.
! Occasionally does a minumum of IO by INQUIREing existance of a
! file.
!
! INTENT (IN)
! character (len=*) :: fit_file, scf_file, ham_file, eig_file
!
implicit none
logical, intent (out) :: recover
integer, intent (out) :: recover_mode
character (len=NCHAR), intent (out) :: rec_file
!** End of interface ***************************************
character (len=13) :: rec_option
recover = .false.
recover_mode = options_recover()
select case (recover_mode)
case (recover_fitcoeff)
recover = .true.
rec_file = fit_file
rec_option = 'read_fitcoeff'
case (recover_scfstate)
recover = .true.
rec_file = scf_file
rec_option = 'read_scfstate'
case (recover_ksmatrix)
recover = .true.
rec_file = ham_file
rec_option = 'read_ksmatrix'
case (recover_eigenvec)
recover = .true.
rec_file = eig_file
rec_option = 'read_eigenvec'
case (recover_fragment)
recover = .true.
rec_file = frg_file
rec_option = 'read_fragment'
end select
if (recover) then
inquire (file=recfile (rec_file), exist=recover)
if (.not. recover) then
if (operations_geo_opt .or. operations_qm_mm) then
recover_mode = recover_nothing
call say ('Warning: No Recover File "' // rec_file // '"')
call say (' "'//trim (rec_option)//'" temporarily turned off !')
else
call error_handler &
('Recover File "' // rec_file // '" not found !')
endif
endif