2
c Main spin-orbit DFT driver
4
logical function dft_scf_so
5
& (rtdb, Etold, Enuc, iVcoul_opt, iVxc_opt,
6
& iter, g_dens, g_dens_at, g_movecs, g_vxc,
7
& g_fock, g_svecs, svals, g_xcinv, g_s)
9
c $Id: dft_scf_so.F 21279 2011-10-24 03:13:09Z niri $
14
integer rtdb ! [input]
15
double precision Etold, Enuc, trace
18
integer iter, swap(20),nswap, ndet, idet
19
integer g_dens(2), g_movecs(2), g_vxc(4),
23
double precision svals(*)
25
integer g_densso(2), g_tmp_ri,
26
& g_moso(2), g_old(2),
27
& g_fockso(2), g_scr2, g_damp_so(2), g_gmovecs(2)
29
integer la, ia ! complex*16 a(nbf_ao, nbf_ao)
30
integer lw, iw ! double precision w(nbf_ao)
32
integer lwork, iwork ! complex*16 work(3)
33
integer lrwork, irwork ! double precision rwork
40
c declarations for fractional occupation
42
integer nmo_fon ! number of fractionally occupied orbitals
43
integer ncore_fon ! number of fully occupied orbitals
44
double precision nel_fon ! fractional electron number
45
double precision avg_fon ! fractional occupancy (averaged)
46
integer nTotOcc ! nTotOcc: no. of (occupied maybe f.o.) mo's
49
double precision rtmp_fon(2), pstrace
50
double precision scale, ncheck
51
integer kfon_occ, lfon_occ
52
logical debug_fon, det_eng
56
double precision rho_n, toll_s
59
double precision ener_scal, ener_kin
60
double precision numelecs
63
character*7 vecs_or_dens
69
logical xc_chkdispauto
70
external xc_chkdispauto
75
#include "mafdecls.fh"
82
#include "zora.fh" ! zora contribution
83
#include "case.fh" ! coulomb attenuation
85
Logical movecs_write_so, movecs_converged
86
External movecs_write_so, movecs_converged
88
Logical movecs_read_header_so, movecs_read_so
89
External movecs_read_header_so, movecs_read_so
93
integer ga_create_atom_blocked
94
external ga_create_atom_blocked
96
logical oprint, status
97
double precision Exc(2), rms(2), derr(2)
98
integer nmo(2), icall(2)
100
integer n3c_dbl, n3c_int, n_batch
102
integer l_3cwhat, k_3cwhat, l_3cERI, k_3cERI
103
integer dft_n3cint, n_semi_bufs, fd
105
double precision dft_n3cdbl
109
integer natoms, nTotEl
112
integer i, j, i1, jstart
114
integer g_tmp, g_fockt, g_wght, g_xyz,g_nq
118
integer nheap, nstack
122
integer itol_max, iaoacc_max
123
integer itol_min, iAOacc_min
124
double precision tol_rho_min, tol_rho_max
126
integer leneval, lcd_coef, icd_coef
127
integer lcntoce, icntoce, lcntobfr, icntobfr,
128
& lcetobfr, icetobfr, lrdens_atom, irdens_atom,
130
double precision start_wall, current_wall, elapsed_wall,
131
& save_wall, current_cpu, start_cpu,
133
integer int_wall_time_reqd
134
double precision ecoul, ecore, noso
135
double precision pp, delta
136
double precision anucl_charg, anel
137
double precision anoca, anocb, onempp
138
double precision etnew, tol2e, tol2e_sleazy
139
c convergence declarations
140
double precision rlshift_input, rlshift_def
141
integer ndamp_input, ndamp_def
143
c Note, damping, levelshifting, and diising logicals
144
c are used to turn on/off these procedures per
145
c iteration. The alternative logicals nodamping,
146
c nolevelshifting, and nodiis are specified and held
147
c for the entire convergence sequence.
149
logical diising, damping, levelshifting
150
logical keep_damp_on,keep_levl_on, keep_diis_on
151
Logical IOLGC, mulliken
152
logical converged, wght_GA
154
logical oprint_parm, oprint_conv, oprint_vecs,
155
& oprint_eval, oprint_syma, oprint_time,
156
& oprint_info, oprint_tol, oprint_final_vecs,
157
& oprint_energy_step, oprint_intermediate_fock,
158
& oprint_3c2e, oprint_interm_overlap, oprint_interm_S2,
159
& oprint_conv_details
160
double precision zero, onem, one, mone
161
parameter(zero = 0.d0, one = 1.d0, mone=-1.0d0, onem = -one)
163
integer ilo, ihi ! For printing movecs analysis
164
double precision eval_pr_tol_lo, eval_pr_tol_hi
165
parameter (eval_pr_tol_lo = -1.5d0, eval_pr_tol_hi=0.5)
168
c early convergence tolerances
170
parameter(itol_min = 7, iAOacc_min = 12, tol_rho_min = 1.d-7)
172
double precision dft_dencvg, dft_time
174
double precision homo, lumo, homo_lumo_gap
176
logical last_time_energy
177
logical check_shift, lmaxov_sv
180
character*255 basis_name, basis_trans
181
integer nopen, nclosed
183
logical cphf_poliz, do_poliz
186
character*255 title1 ! Returns title of job that created vectors
187
character*255 basis_name1 ! Returns name of basis set
188
character*255 scftype1 ! Returns the SCF type of the vectors
189
integer nbf1 ! Returns no. of functions in basis
194
logical dft_zora_read_so, dft_zora_write_so
195
external dft_zora_read_so, dft_zora_write_so
197
logical dft_zora_inquire_file_so
198
external dft_zora_inquire_file_so
200
character*255 zorafilename
202
integer g_zora_scale_sf(2)
204
integer g_zora_scale_so(3)
205
double precision Ezora_sf
206
integer switch_sclMO_so ! switch =1,0 ON,OFF sclMO-FA-02-18-11
211
external dft_scaleMO_so ! FA-occupations from input script 02-15-11
219
c to do scalar relativistic calculations within the two-component framework
221
do_purescalar = .false.
222
if (.not.rtdb_get(rtdb,'sodft:scalar',mt_log,1,do_purescalar))
223
& do_purescalar = .false.
225
if (me.eq.0.and.do_purescalar) then
226
call util_print_centered(LuOut,
227
$ 'Neglecting spin-orbit terms', 23, .true.)
231
call ecce_print_module_entry('dft')
235
oprint = util_print('information', print_low)
236
oprint_info = util_print('common', print_debug)
237
oprint_parm = util_print('parameters', print_default)
238
oprint_3c2e = util_print('3c 2e integrals', print_default)
239
oprint_conv = util_print('convergence', print_default)
240
oprint_conv_details = util_print('convergence details',
242
oprint_vecs = util_print('intermediate vectors', print_high)
243
oprint_eval = util_print('intermediate evals', print_high)
244
oprint_syma = util_print('interm vector symm', print_high)
245
oprint_time = util_print('dft timings', print_high)
246
oprint_tol = util_print('screening parameters', print_high)
247
oprint_energy_step = util_print('intermediate energy info',
249
oprint_intermediate_fock = util_print('intermediate fock matrix',
251
oprint_interm_S2 = util_print('intermediate S2',print_high)
252
oprint_interm_overlap = util_print('intermediate overlap',
254
oprint_final_vecs = util_print('final vectors', print_high)
257
call int_1e_uncache_ga()
259
c Store SCF hamiltonian type as DFT for use in BGJ routines
260
if (.not. rtdb_put(rtdb, 'bgj:scf_type', MT_INT, 1, 2))
261
$ call errquit('dft_scf_so: put of bgj:scf_type failed',0,
265
c see if levelshifting monitoring is desired
267
if (.not. rtdb_get(rtdb, 'dft:check_shift', mt_log, 1,
269
check_shift = .false.
272
if (.not. geom_ncent(geom, natoms))
273
& call errquit('dft_scf_so: geom_ncent failed',73, GEOM_ERR)
274
if (.not. geom_nuc_charge(geom, anucl_charg))
275
& call errquit('dft_scf_so: geom_nuc_charge failed',
278
anel = int(anucl_charg) - rcharge
280
c Pre-compute mapping vectors
283
& (mt_int,nshells_ao,'cntoce map',lcntoce,icntoce))
284
& call errquit('dft_scf_so:push_get failed', 13, MA_ERR)
286
& (mt_int,nshells_ao*2,'cntoce map',lcntobfr,icntobfr))
287
& call errquit('dft_scf_so:push_get failed', 13, MA_ERR)
289
& (mt_int,natoms*2,'cntoce map',lcetobfr,icetobfr))
290
& call errquit('dft_scf_so:push_get failed', 13, MA_ERR)
292
call build_maps(ao_bas_han, int_mb(icntoce), int_mb(icntobfr),
293
& int_mb(icetobfr), natoms, nshells_ao)
295
c Set aside some memory for reduced density matrix
297
if (.not.MA_Push_Get(MT_Dbl,2*natoms*natoms,'rdens_atom',
298
& lrdens_atom,irdens_atom))
299
& call errquit('dft_scf_so: cannot allocate rdens_atom',
302
c determine pattern of orbitals' occupancy
304
if (.not. MA_Push_Get(MT_Dbl,nbf_ao*2,'mo occ',l_occ,k_occ))
305
& call errquit('dft_scf_so: failed to alloc',999, MA_ERR)
307
c get orbital overlap tolerance
309
if (.not. rtdb_get(rtdb, 'dft:toll_s', MT_DBL, 1, toll_s))
310
. call errquit('dft_scf_so: lost toll_s ',0, RTDB_ERR)
312
nTotEl = noc(1) + noc(2)
321
call dfill(nbf_mo, 0.0d0, dbl_mb(k_occ), 1)
323
dbl_mb(i-1+k_occ) = 1.0d0
325
do i = nbf_ao+1, nbf_ao+noc(2)
326
dbl_mb(i-1+k_occ) = 1.0d0
331
c Determine whether to fit the electronic charge density.
334
if (iVcoul_opt.eq.1)CDFIT = .TRUE.
336
if (iVxc_opt.eq.1)XCFIT = .TRUE.
338
c Define various constants.
340
npol = (ipol*(ipol+1))/2
344
tol_rho_max = tol_rho
346
& call dft_tstamp(' Before 3c-2e initialize.')
350
I natoms,npol,oprint_parm,oprint_3c2e,
351
O n3c_int,n3c_dbl,n_semi_bufs,
352
O l_3ceri,k_3ceri, l_3cwhat,k_3cwhat)) then
353
call dft_3cincor(n_batch, n3c_int, int_mb(k_3cwhat),
354
& dbl_mb(k_3cERI), n3c_dbl, iwhat_max,
358
if (me.eq.0 .and. oprint_3c2e)write(LuOut,3230)
362
3230 format(/,10x,'Incore memory use for 3-center 2e- integrals is ',
365
if (imull.eq.1)mulliken = .true.
367
if (noio.eq.1)IOLGC = .FALSE.
369
c Energy decomposition switch
374
c SCF energy convergence criterion.
376
g_tmp = ga_create_atom_blocked(geom, AO_bas_han, 'ga_temp')
377
g_fockt = ga_create_atom_blocked(geom, AO_bas_han, 'fock tr')
379
c Set up local convergence parameters
383
levelshifting = levelshift
384
keep_damp_on = .false.
385
keep_levl_on = .false.
386
keep_diis_on = .false.
388
rlshift_input = rlshift
391
rlshift = rlshift_def
393
if (nodamping)damping = .false.
394
if (nolevelshifting) then
395
levelshifting = .false.
396
rlshift = rlshift_def
406
levelshifting = .true.
407
rlshift = rlshift_input
413
c Initialize DIIS call counter.
418
c Begin the SCF cycle.
420
c allocate eigenvalue array
423
if (.not.MA_Push_Get(MT_Dbl,leneval,'eval',l_eval,k_eval(1)))
424
& call errquit('dft_scf_so: cannot allocate eval',0, MA_ERR)
425
k_eval(2) = k_eval(1) + nbf_mo
427
c Dump DFT parameters (if debugging) to see if they make sense
429
if (me.eq.0.and.oprint_info)call dft_dump_info(me)
431
c Get initial density.
434
& call dft_tstamp(' Before call to DFT_INIT.')
437
c allocate array for irreps
439
if (.not.MA_Push_Get(mt_int,2*nbf_ao,'dft:irreps',l_ir,k_ir))
440
& call errquit('dft_scf_so: cannot allocate irreps',0, MA_ERR)
442
nclosed = (nTotEl - nopen) / 2
444
if (.not. bas_name(ao_bas_han, basis_name, basis_trans))
445
$ call errquit('dft_scf_so: bas_name?', 0, BASIS_ERR)
447
c get info for int2e_ and set sleazy tolerance
450
call scf_get_fock_param(rtdb, tol2e_sleazy)
452
c Force sleazy SCF into "direct" mode.
454
call fock_force_direct(rtdb)
456
cso allocate Fock matrix and movecs
459
c real molecular orbital vectors
460
if(.not.ga_create(mt_dbl,nbf_mo,nbf_mo,'Movecs Re',0,0,
462
& call errquit('dft_scf_so: error creating Movecs Re',0,
464
if(.not.ga_create(mt_dbl,nbf_ao,nbf_ao,'Movecs Re',0,0,
466
& call errquit('dft_scf_so: error creating Movecs Re',0,
469
c imaginary molecular orbital vectors
470
if(.not.ga_create(mt_dbl,nbf_mo, nbf_mo,'Movecs Im',0,0,
472
& call errquit('dft_scf_so: error creating Movecs Im',0,
474
if(.not.ga_create(mt_dbl,nbf_ao, nbf_ao,'Movecs Im',0,0,
476
& call errquit('dft_scf_so: error creating Movecs Im',0,
479
call ga_zero(g_moso(1))
480
call ga_zero(g_moso(2))
481
call ga_zero(g_gmovecs(1))
482
call ga_zero(g_gmovecs(2))
486
c real part of the fock matrix
487
if(.not.ga_create(mt_dbl,nbf_mo,nbf_mo,'Fock Re',0,0,
489
& call errquit('dft_scf_so: error creating Fock Re',0, GA_ERR)
490
call ga_zero(g_fockso(1))
492
c imaginary part of the fock matrix
493
if(.not.ga_create(mt_dbl,nbf_mo,nbf_mo,'Fock Im',0,0,
495
& call errquit('dft_scf_so: error creating Fock Im',0, GA_ERR)
496
call ga_zero(g_fockso(2))
499
if(.not.ga_create(mt_dbl,nbf_mo,nbf_mo,'old re',0,0,
501
& call errquit('dft_scf_so: error creating Old Re',0, GA_ERR)
502
call ga_zero(g_old(1))
503
if(.not.ga_create(mt_dbl,nbf_mo,nbf_mo,'old Im',0,0,
505
& call errquit('dft_scf_so: error creating Old Im',0, GA_ERR)
506
call ga_zero(g_old(2))
508
c real part of the spin-orbit density matrix
509
if(.not.ga_create(mt_dbl,nbf_mo,nbf_mo,'DenMx Re',0,0,
511
& call errquit('dft_scf_so: error creating DenMx Re',0, GA_ERR)
512
call ga_zero(g_densso(1))
514
c imaginary part of the spin-orbit density matrix
515
if(.not.ga_create(mt_dbl,nbf_mo,nbf_mo,'DenMx Im',0,0,
517
& call errquit('dft_scf_so: error creating DenMx Im',0, GA_ERR)
518
call ga_zero(g_densso(2))
521
if(.not.ga_create(mt_dbl,nbf_mo,nbf_mo,'Tmp ReIm',0,0,
523
& call errquit('dft_scf_so: error creating Tmp ReIm',0, GA_ERR)
524
call ga_zero(g_tmp_ri)
526
c spin-orbit matrices: 1->z, 2->y, 3->x
527
if(.not.ga_create(mt_dbl,nbf_ao,nbf_ao,'so z',0,0,
529
& call errquit('dft_scf_so: error creating so z',0, GA_ERR)
530
call ga_zero(g_so(1))
531
if(.not.ga_create(mt_dbl,nbf_ao,nbf_ao,'so y',0,0,
533
& call errquit('dft_scf_so: error creating so y',0, GA_ERR)
534
call ga_zero(g_so(2))
535
if(.not.ga_create(mt_dbl,nbf_ao,nbf_ao,'so x',0,0,
537
& call errquit('dft_scf_so: error creating so x',0, GA_ERR)
538
call ga_zero(g_so(3))
541
if(.not.ga_create(mt_dbl, 2*nbf, 2*nbf,'old den', 0, 0,
543
& call errquit('dft_scf_so: error creating damp ga', 0, GA_ERR)
544
call ga_zero(g_damp_so(1))
545
if(.not.ga_create(mt_dbl, 2*nbf, 2*nbf,'old den', 0, 0,
547
& call errquit('dft_scf_so: error creating damp ga', 0, GA_ERR)
548
call ga_zero(g_damp_so(2))
553
call util_print_centered(LuOut,
554
$ 'Performing ZORA calculations', 23, .true.)
558
c == get filename for zora data ==
559
call util_file_name('zora_so',.false.,.false.,zorafilename)
561
c == zora: scalar arrays ==
562
if(.not.ga_create(mt_dbl,nbf_ao,nbf_ao,'g_zora_sf',0,0,
564
& call errquit('dft_scf_so: error creating g_zora_sf',0,
566
call ga_zero(g_zora_sf(1))
567
if(.not.ga_create(mt_dbl,nbf_ao,nbf_ao,'g_zora_sf',0,0,
569
& call errquit('dft_scf_so: error creating g_zora_sf',0,
571
call ga_zero(g_zora_sf(2))
573
c == zora: scalar energy scaling arrays ==
574
if(.not.ga_create(mt_dbl,nbf_ao,nbf_ao,'g_zora_scale_sf',0,0,
575
& g_zora_scale_sf(1)))
576
& call errquit('dft_scf_so: error creating g_zora_scale_sf',0,
578
call ga_zero(g_zora_scale_sf(1))
579
if(.not.ga_create(mt_dbl,nbf_ao,nbf_ao,'g_zora_scale_sf',0,0,
580
& g_zora_scale_sf(2)))
581
& call errquit('dft_scf_so: error creating g_zora_scale_sf',0,
583
call ga_zero(g_zora_scale_sf(2))
585
c == zora: spin-orbit arrays ==
586
if(.not.ga_create(mt_dbl,nbf_ao,nbf_ao,'g_zora_so 1',0,0,
588
& call errquit('dft_scf_so: error creating g_zora_so 1',0, GA_ERR)
589
call ga_zero(g_zora_so(1))
590
if(.not.ga_create(mt_dbl,nbf_ao,nbf_ao,'g_zora_so 2',0,0,
592
& call errquit('dft_scf_so: error creating g_zora_so 2',0, GA_ERR)
593
call ga_zero(g_zora_so(2))
594
if(.not.ga_create(mt_dbl,nbf_ao,nbf_ao,'g_zora_so 3',0,0,
596
& call errquit('dft_scf_so: error creating g_zora_so 3',0, GA_ERR)
597
call ga_zero(g_zora_so(3))
599
c == zora: spin-orbit energy scaling arrays ==
600
if(.not.ga_create(mt_dbl,nbf_ao,nbf_ao,'g_zora_scale_so 1',0,0,
601
& g_zora_scale_so(1)))
602
& call errquit('dft_scf_so: error creating g_zora_scale_so 1',0,
604
call ga_zero(g_zora_scale_so(1))
605
if(.not.ga_create(mt_dbl,nbf_ao,nbf_ao,'g_zora_scale_so 2',0,0,
606
& g_zora_scale_so(2)))
607
& call errquit('dft_scf_so: error creating g_zora_scale_so 2',0,
609
call ga_zero(g_zora_scale_so(2))
610
if(.not.ga_create(mt_dbl,nbf_ao,nbf_ao,'g_zora_scale_so 3',0,0,
611
& g_zora_scale_so(3)))
612
& call errquit('dft_scf_so: error creating g_zora_scale_so',0,
614
call ga_zero(g_zora_scale_so(3))
616
c == generate an superposition of atomic densities ==
617
call ga_zero(g_dens_at(1))
618
if (ipol.gt.1) call ga_zero(g_dens_at(2))
619
call guess_dens(rtdb, geom, ao_bas_han, g_dens_at)
620
if (oskel) call ga_symmetrize(g_dens_at(1))
622
call ga_copy(g_dens_at(1),g_dens_at(2))
623
call ga_dscal(g_dens_at(1),dble(ntotel-nclosed)/(ntotel))
624
call ga_dscal(g_dens_at(2),dble(nclosed)/(ntotel))
625
if(oskel) call ga_symmetrize(g_dens_at(2))
628
c == in case fon is used together with zora ==
629
c == pstrace is queried in the grid code ==
632
if (rtdb_get(rtdb,'dft:fon',mt_log,1,fon)) then
633
pstrace=ga_ddot(g_dens_at(1),g_s)
634
pstrace=pstrace + ga_ddot(g_dens_at(2),g_s)
635
if(ga_nodeid().eq.0) write (luout,'(5x,a,1x,e15.7)')
636
& 'tr(P*S): ',pstrace
637
if (.not. rtdb_put(rtdb, 'dft:pstrace', mt_dbl, 1, pstrace))
638
& call errquit('dft_scf: rtdb_put pstrace failed', 1, RTDB_ERR)
641
c == try reading the zora contributions from file ==
643
if (.not.dft_zora_read_so(zorafilename, nbf_ao, ipol, nmo, mult,
644
& g_zora_sf, g_zora_scale_sf, g_zora_so, g_zora_scale_so))
647
if (icalczora.eq.1) then
649
call util_print_centered(LuOut,
650
$ 'Generating atomic ZORA corrections', 23, .true.)
654
c == calculate the zora atomic corrections ==
655
call zora_getv_so(rtdb, g_dens_at, g_zora_sf, g_zora_scale_sf,
656
& g_zora_so, g_zora_scale_so, nexc)
658
c == write out the atomic zora corrections to file ==
659
if (.not.dft_zora_write_so(rtdb, ao_bas_han, zorafilename,
660
& nbf_ao, ipol, nmo, mult, g_zora_sf, g_zora_scale_sf,
661
& g_zora_so, g_zora_scale_so))
662
& call errquit('dft_scf_so: dft_zora_write_so failed', 0,
667
c == for scalar calculations via the 2 component formalism
668
if (do_purescalar) then
669
call ga_zero(g_zora_so(1))
670
call ga_zero(g_zora_so(2))
671
call ga_zero(g_zora_so(3))
672
call ga_zero(g_zora_scale_so(1))
673
call ga_zero(g_zora_scale_so(2))
674
call ga_zero(g_zora_scale_so(3))
679
c check for fon input
682
if (rtdb_get(rtdb,'dft:fon',mt_log,1,fon)) then
684
c variable 'fon' should be true, otherwise there's
685
c something fishy going on:
686
if (.not.fon) call errquit(
687
& 'dft_scf_so: fon stored in RTDB but not .true.', 1,
690
c note: *_fon variables are read here not as arrays with
691
c two elements, assuming that the input didn't specify
692
c spects for alpha and beta separately
694
if (.not.rtdb_get(rtdb,'dft:nmo_fon',mt_int,2,ntmp_fon)) then
696
write(LuOut,*)"Error: fractional occupation
697
& calculation specified without setting number of orbitals"
699
call errquit('dft_scf_so: nmo_fon rtdb_get failed', 1,
702
nmo_fon = ntmp_fon(1)
704
if (.not.rtdb_get(rtdb,'dft:nel_fon',mt_dbl,2,rtmp_fon)) then
706
write(LuOut,*)"Error: fractional occupation
707
& specified without setting number of electrons"
709
call errquit('dft_scf_so: nel_fon rtdb_get failed', 1,
712
nel_fon = rtmp_fon(1)
714
if (.not.rtdb_get(rtdb,'dft:ncore_fon',mt_int,2,ntmp_fon)) then
716
write(LuOut,*)"Error: fractional occupation
717
& calculation specified without setting filled levels"
719
call errquit('dft_scf_so: nel_core rtdb_get failed', 1,
722
ncore_fon = ntmp_fon(1)
724
if (rtdb_get(rtdb, 'dft:debugfon', mt_log, 1,
725
& debug_fon)) continue
727
else ! keyword not in RTDB
733
call ga_zero(g_densso(1))
734
call ga_zero(g_densso(2))
736
spinor=spinor_guess(movecs_in)
740
c fractional occupations:
742
c ... jochen: presumably good enough for initial guess. We won't tinker
743
c with that. later, the fractional occupations
744
c are caluclated explicitly
747
noc(1)=(noc(1)+noc(2)-nel_fon)/2
748
noc(2)=(noc(1)+noc(2)-nel_fon)/2
751
call dft_guessin(movecs_in,ldmix,ncanorg,fon,
752
& vecs_or_dens, ipol,nbf_ao,g_movecs,g_gmovecs,
755
call scf_vectors_guess(rtdb, tol2e_sleazy, geom, ao_bas_han,
756
& basis_trans, movecs_in, movecs_out,
757
& movecs_guess, scftype, nclosed, nopen,
758
& nbf, nmo, noc(1), noc(2), k_eval, k_occ,
759
& k_ir, g_gmovecs, g_dens, 'density',
760
& 'dft', title, oskel, oadapt,
763
call dft_guessout(nmo,nbf_ao,g_gmovecs,g_movecs,ipol)
765
c fon: undo temp setting of noc(:)
773
call dfill(nbf_mo, 0.0d0, dbl_mb(k_occ), 1)
775
dbl_mb(i-1+k_occ) = 1.0d0
778
c map initial guess movecs from spin-free calculations g_moso(1)
779
c noc(1).ge.noc(2) is assumed
781
do i=1,min(noc(1),noc(2))
782
call ga_dadd_patch(1.d0,g_movecs(1),1,nbf_ao,i,i,
783
$ 0.d0,g_moso(1),1,nbf_ao,2*(i-1)+1,2*(i-1)+1,
784
$ g_moso(1),1,nbf_ao,2*(i-1)+1,2*(i-1)+1)
785
call ga_dadd_patch(1.d0,g_movecs(2),1,nbf_ao,i,i,
786
$ 0.d0,g_moso(1),1+nbf_ao,nbf_mo,2*(i-1)+2,2*(i-1)+2,
787
$ g_moso(1),1+nbf_ao,nbf_mo,2*(i-1)+2,2*(i-1)+2)
790
call ga_dadd_patch(1.d0,g_movecs(1),1,nbf_ao,i,i,
791
$ 0.d0,g_moso(1),1,nbf_ao,noc(2)+i,noc(2)+i,
792
$ g_moso(1),1,nbf_ao,noc(2)+i,noc(2)+i)
795
call ga_dadd_patch(1.d0,g_movecs(2),1,nbf_ao,i,i,
796
$ 0.d0,g_moso(1),1+nbf_ao,nbf_mo,noc(1)+i,noc(1)+i,
797
$ g_moso(1),1+nbf_ao,nbf_mo,noc(1)+i,noc(1)+i)
800
call ga_dadd_patch(1.d0,g_movecs(1),1,nbf_ao,i,i,
801
$ 0.d0,g_moso(1),1,nbf_ao,2*(i-1)+1,2*(i-1)+1,
802
$ g_moso(1),1,nbf_ao,2*(i-1)+1,2*(i-1)+1)
803
call ga_dadd_patch(1.d0,g_movecs(2),1,nbf_ao,i,i,
804
$ 0.d0,g_moso(1),1+nbf_ao,nbf_mo,2*(i-1)+2,2*(i-1)+2,
805
$ g_moso(1),1+nbf_ao,nbf_mo,2*(i-1)+2,2*(i-1)+2)
811
c read spinors from files
813
c get MO vectors from file
815
c if (.not. rtdb_cget(rtdb, 'dft:input vectors', 1, movecs_in))
816
c $ call errquit('dft_scf_so: DFT MO vectors not defined',0)
817
status = movecs_read_header_so(movecs_in, title1, basis_name1,
820
c Should check much more info than just nbf for consistency
825
if (2*nbf_ao .ne. nbf1)then
826
write(6,*)'dft_scf_so movecs output = ',movecs_in
827
call errquit('dft_scf_so: could not read mo vectors',911,
831
status = status .and.
832
$ movecs_read_so(movecs_in, dbl_mb(k_occ),
833
$ dbl_mb(k_eval(1)), g_moso)
837
write(6,*)'dft_scf_so movecs output = ',movecs_in
838
call errquit('dft_scf_so: could not read mo vectors',917,
842
call movecs_swap_so(rtdb,'dft',scftype,g_moso,
843
& dbl_mb(k_occ),dbl_mb(k_eval(1)))
846
c Form Re and Im of density matrix
848
c ... jochen 10/11: implemented new way of applying
849
c fractional occupations. See dft_densm.F of the scalar
852
c$$$ if (.not.rtdb_get(rtdb,'sodft:fon',mt_log,1,fon))
854
c$$$ if (.not.rtdb_get(rtdb,'sodft:nmo_fon',mt_int,1,nmo_fon))
856
c$$$ if (.not.rtdb_get(rtdb,'sodft:nel_fon',mt_int,1,nel_fon))
858
c$$$ nTotOcc = (nTotEl-nel_fon) + nmo_fon
860
nTotOcc = nTotEl ! if there is no fractional occupation
862
c$$$ avg_fon = dble(nel_fon)/dble(nmo_fon)
863
c$$$ do i = (nTotEl-nel_fon)+1, nTotOcc
864
c$$$ dbl_mb(i-1+k_occ) = avg_fon
867
if (nmo_fon.lt.1) call errquit(
868
& 'dft_scf_so:fon nmo_fon <1',
870
if (nel_fon.lt.0d0) call errquit(
871
& 'dft_scf_so:fon nel_fon <0',
874
avg_fon = nel_fon/dble(nmo_fon)
875
nTotOcc = ncore_fon + nmo_fon
878
c$$$ write (luout,*) 'DEBUG: occupations before applying FON'
880
c$$$ write (luout,*) i, dbl_mb(i-1+k_occ), dbl_mb(i-1+nbf_ao+k_occ)
885
if (i> 2*nbf_ao) call errquit(
886
& 'dft_densm:fon focc index exceeds 2nbf error 1',
888
dbl_mb(i-1+k_occ) = 1d0
889
ncheck = ncheck + 1d0
891
do i = ncore_fon + 1, ncore_fon + nmo_fon
892
if (i> 2*nbf_ao) call errquit(
893
& 'dft_densm:fon focc index exceeds 2nbf error 2',
895
dbl_mb(i-1+k_occ) = avg_fon
896
ncheck = ncheck + avg_fon
899
if(abs(ncheck-dble(nTotEl)).gt.1d-3 .and. me.eq.0) then
900
write(luout,*) ' frac. electrons ',ncheck,' vs ',nTotEl
903
c$$$ write (luout,*) 'DEBUG: occupations AFTER applying FON'
905
c$$$ write (luout,*) i, dbl_mb(i-1+k_occ), dbl_mb(i-1+nbf_ao+k_occ)
909
switch_sclMO_so=0 ! FA-09-26-11
911
c the fractionally occupied mo's are scaled by the sqrt of the fon
914
do i = ncore_fon + 1, ncore_fon + nmo_fon
915
if (i> nbf_mo) call errquit(
916
& 'dft_densm:fon g_moso index exceeds nbf_mo',
918
scale = sqrt(dbl_mb(i-1+k_occ)) ! sqrt(occupation number)
919
call ga_scale_patch(g_moso(1),
920
& 1, nbf_mo, i, i, scale)
921
call ga_scale_patch(g_moso(2),
922
& 1, nbf_mo, i, i, scale)
924
if(me.eq.0) write(luout,'(5x,a)') 'FON applied'
927
c ---- FA-02-18-11 : occupations keyword ---- START
928
c call dft_scaleMO_so(rtdb,g_moso,dbl_mb(k_occ),g_densso,
929
c & nbf_mo,nTotOcc,switch_sclMO_so)
932
c calculate the spin-orbit density matrix using the scaled mo's
934
c if (switch_sclMO_so.ne.1) then
935
c if (ga_nodeid().eq.0)
936
c & write(*,*) 'ENTER calc so-density matrix std-way'
937
call dft_densm_so(g_densso, g_moso, nbf_ao, nTotOcc)
940
c restore the scaled mo's
943
do i = ncore_fon + 1, ncore_fon + nmo_fon
944
if (i> nbf_mo) call errquit(
945
& 'dft_densm:fon g_moso index exceeds nbf_mo',
947
if (dbl_mb(i-1+k_occ) < 1d-4) call errquit(
948
& 'dft_densm:fon frac occup < 1E-4. Aborting suspisciously',
950
scale = 1d0/sqrt(dbl_mb(i-1+k_occ)) ! 1/sqrt(occupation number)
951
call ga_scale_patch(g_moso(1),
952
& 1, nbf_mo, i, i, scale)
953
call ga_scale_patch(g_moso(2),
954
& 1, nbf_mo, i, i, scale)
958
c calculate the spin-free density matrix from the spin-orbit density matrix
960
call ga_zero(g_dens(1))
961
call ga_zero(g_dens(2))
962
call ga_dens_sf(g_dens, g_densso(1), nbf_ao)
964
c check fon, calculate tr[P S] and print
965
c if (debug_fon) call dft_pstrace(g_dens,ao_bas_han,nbf_ao,oskel)
967
pstrace=ga_ddot(g_dens(1),g_s)
968
pstrace=pstrace + ga_ddot(g_dens(2),g_s)
969
if(ga_nodeid().eq.0) write (luout,'(5x,a,1x,e15.7)')
970
& 'tr(P*S): ',pstrace
971
if (.not. rtdb_put(rtdb, 'dft:pstrace', mt_dbl, 1, pstrace))
972
& call errquit('dft_scf: rtdb_put pstrace failed', 1, RTDB_ERR)
977
call fock_2e_tidy(rtdb)
979
c set initial coulomb acc
981
c write(6,*)' movecs_guess = ',movecs_guess
982
if (movecs_guess.eq.'restart')ltight=.true.
984
c May not want levelshifting initially until sure that the
985
c transformed Fock matrix will be diagonally dominant, or
986
c alternatively shift the piss out of it.
988
if (movecs_guess.eq.'restart')then
989
levelshifting = .true.
991
levelshifting = .false.
999
tol_rho = tol_rho_max
1002
itol2e = min(itol_min,itol_max)
1003
iAOacc = min(iAOacc_min,iAOacc_max)
1004
tol_rho = max(tol_rho_min,tol_rho_max)
1007
tol2e = 10.d0**(-itol_max)
1009
c Restore SCF parameters
1011
call scf_get_fock_param(rtdb, tol2e)
1013
c If open shell put the total density matrix in g_dens(1)
1015
call ga_dadd(one,g_dens(1),one,g_dens(2),g_dens(1))
1018
c Call to Mulliken Pop Analysis for initial density
1021
if (me.eq.0)call dft_header
1022
& (' Total Density - Mulliken Population Analysis')
1023
call mull_pop(geom,ao_bas_han,g_dens(1),g_s,'total')
1025
c analysis of spin density
1027
if (me.eq.0) call dft_header
1028
& (' Spin Density - Mulliken Population Analysis')
1029
call ga_dadd(one,g_dens(1),-2.d0,g_dens(2),g_dens(2))
1030
call mull_pop(geom,ao_bas_han,g_dens(2),g_s,'spin')
1032
c restore beta density in g_dens(2)
1034
call ga_dadd(one,g_dens(1),-1.d0,g_dens(2),g_dens(2))
1035
call ga_dscal(g_dens(2),0.5d0)
1040
cso ma for complex diagonalizer
1041
if (.not.MA_Push_Get(MT_DCpl,nbf_mo*nbf_mo,'cpl a',la,ia))
1042
& call errquit('dft_scf: cannot allocate cpl a',0, MA_ERR)
1043
if (.not.MA_Push_Get(MT_Dbl,nbf_mo,'cpl eval',lw,iw))
1044
& call errquit('dft_scf: cannot allocate cpl eval',0, MA_ERR)
1045
llwork = max(1, 2*nbf_mo-1)
1046
if (.not.MA_Push_Get(MT_DCpl,llwork,'cpl work',lwork,iwork))
1047
& call errquit('dft_scf: cannot allocate cpl work',0, MA_ERR)
1048
if (.not.MA_Push_Get(MT_Dbl,max(1,3*nbf_mo-2),'w.s',lrwork,
1050
& call errquit('dft_scf: cannot allocate w.s',0, MA_ERR)
1051
if (.not.ma_push_get(mt_dbl,nbf_mo,'buff',lbuff,ibuff))
1052
& call errquit('dft_scf:push_get failed', 13, MA_ERR)
1056
c Top of infinite SCF iteration loop
1058
c Write prep time required
1061
if (me.eq.0.and.oprint)then
1062
current_cpu = util_cpusec()
1063
write(LuOut,20)current_cpu
1064
20 format(2x,' Time prior to 1st pass: ',f8.1)
1067
c start DFT_SCF timer
1069
start_wall = util_wallsec()
1070
start_cpu = util_cpusec()
1071
dft_time = -start_cpu
1074
& call dft_tstamp(' Before SCF iter loop. ')
1076
last_time_energy = .false.
1080
1000 continue ! top of the scf loop
1082
if (me.eq.0 .and. oprint_conv_details)
1083
& write(LuOut,124)damping, levelshifting, diising
1084
124 format(10x,' DAMPING=',l1,' LEVELSHIFTING=',l1,
1087
if (me.eq.0.and.oprint_tol)write(LuOut,3234)itol2e,iAOacc,iXCacc
1088
3234 format(10x,'itol2e=',i2,' iAOacc=',i2,' iXCacc=',i2)
1094
call ga_zero(g_fockso(1))
1095
call ga_zero(g_fockso(2))
1097
c Accumulate core hamiltonian into Fock matrix;
1098
c compute core energy
1100
call ga_zero(g_fock)
1101
call int_1e_ga(ao_bas_han, ao_bas_han, g_fock, 'kinetic', oskel)
1103
call int_1e_ga(ao_bas_han, ao_bas_han, g_fock, 'potential', oskel)
1104
call ga_fock_sf(g_fock, g_fockso(1), nbf_ao)
1106
cso Re(Dsf)=Re(Daa)+Re(Dbb)=g_dens(1)
1107
cso <Hsf> = Re(Dsf) dot T+Vsf
1108
Ecore = ga_ddot(g_dens(1), g_fock)
1110
cso Accumulate s.o. contribution to fock matrix
1113
call ga_zero(g_so(1))
1114
call ga_zero(g_so(2))
1115
call ga_zero(g_so(3))
1117
cso Calculate the spin-orbit contributions from the ecp
1118
if( .not. do_zora .and. .not. do_purescalar) then
1119
call int_1e_ga(ao_bas_han, ao_bas_han, g_so, 'so', oskel)
1120
call ga_scale(g_so(1),dble(0.5d0)) ! z
1121
call ga_scale(g_so(2),dble(0.5d0)) ! y
1122
call ga_scale(g_so(3),dble(0.5d0)) ! x
1124
cso Add in the s.o. contribution to the fock matrix
1125
call ga_fock_so(g_so, g_fockso, nbf_ao)
1128
cso Accumulate z-component s.o. contribution
1129
cso Re(Dz)=-Im(Daa)+Im(Dbb)
1130
cso <Hz>=Re(Dz) dot Vz
1133
call ga_dens_so(g_tmp, g_densso, nbf_ao, 'z')
1134
Ecore = Ecore + ga_ddot(g_tmp, g_so(1))
1135
! write(*,*)"Ecore+=so(z)", ecore
1137
c == add in the spin-orbit zora contribution (z) ==
1138
if (do_zora) Ecore = Ecore + ga_ddot(g_tmp, g_zora_so(1))
1140
cso Accumulate y-component s.o. contribution
1141
cso Re(Dy)=Re(Dab)-Re(Dba)
1142
cso <Hy>=Re(Dy) dot Vy
1145
call ga_dens_so(g_tmp, g_densso, nbf_ao, 'y')
1146
Ecore = Ecore + ga_ddot(g_tmp, g_so(2))
1147
! write(*,*)"Ecore+=so(y)", ecore
1149
c == add in the spin-orbit zora contribution (y) ==
1150
if (do_zora) Ecore = Ecore + ga_ddot(g_tmp, g_zora_so(2))
1152
cso Accumulate x-component s.o. contribution
1153
cso Re(Dx)=-Im(Dab)-Im(Dba)
1154
cso <Hx>=Re(Dx) dot Vx
1157
call ga_dens_so(g_tmp, g_densso, nbf_ao, 'x')
1158
Ecore = Ecore + ga_ddot(g_tmp, g_so(3))
1159
! write(*,*)"Ecore+=so(x)", ecore
1161
c == add in the spin-orbit zora contribution (x) ==
1162
if (do_zora) Ecore = Ecore + ga_ddot(g_tmp, g_zora_so(3))
1166
c Pre-compute reduced total density matrices over atoms
1168
call dfill(ipol*natoms*natoms, 0.0d0, dbl_mb(irdens_atom), 1)
1169
nscr = nbf_ao_mxnbf_ce*nbf_ao_mxnbf_ce
1170
if (.not.MA_Push_Get(MT_Dbl,nscr,'scr',lscr,iscr))
1171
& call errquit('dft_scf: cannot allocate scr',0, MA_ERR)
1172
call util_ga_mat_reduce(nbf_ao, natoms, int_mb(icetobfr), g_dens,
1173
& ipol, dbl_mb(irdens_atom), 'rms',
1174
& dbl_mb(iscr), nbf_ao_mxnbf_ce,.true.)
1175
c write(*,'("irdens",5f10.7)')
1176
c & (dbl_mb(irdens_atom+i),i=0,ipol*natoms*natoms-1)
1177
if (.not.ma_pop_stack(lscr))
1178
& call errquit('dft_scf: cannot pop stack:lscr',0, MA_ERR)
1184
if (cam_exch) call case_setflags(.false.)
1186
c Fit the electron charge density.
1188
if (.not.MA_Push_Get(MT_Dbl,nbf_cd,'cd_coef',lcd_coef,
1190
& call errquit('dft_scf: cannot allocate cd_coef',0, MA_ERR)
1192
& call dft_tstamp(' Before call to FITCD. ')
1193
call dft_fitcd(1,Dbl_MB(icd_coef), dbl_mb(k_3cERI), Ecoul,
1194
& g_dens, nTotEl, n_batch, n3c_int,
1195
& int_mb(k_3cwhat), n3c_dbl, iwhat_max,
1196
& n_semi_bufs, fd, IOLGC,
1198
& .false., 0d0, .false.)
1202
& call dft_tstamp(' Before call to GETVCOUL.')
1203
call dft_getvc(Dbl_MB(icd_coef), dbl_mb(k_3cERI), Ecoul,
1204
& g_tmp, iVcoul_opt, n_batch,
1205
& n3c_int, int_mb(k_3cwhat), n3c_dbl, iwhat_max,
1206
& n_semi_bufs, fd, IOLGC,
1209
c Add V coul to Fock Matrix
1211
cso call ga_dadd(one, g_tmp, one, g_fock, g_fock)
1212
call ga_fock_sf(g_tmp, g_fockso(1), nbf_ao)
1214
if (.not.ma_pop_stack(lcd_coef))
1215
& call errquit('dft_scf: cannot pop stacklcd_coef',0,
1219
c Restore alpha and beta densities.
1221
call ga_dadd(one, g_dens(1), onem, g_dens(2), g_dens(1))
1223
c Note that g_dens(1) now contains the alpha density
1224
c matrix and g_dens(2) contains the beta
1226
c Pre-compute reduced alpha and beta density matrices over atoms
1228
call dfill(ipol*natoms*natoms, 0.0d0, dbl_mb(irdens_atom), 1)
1229
nscr = nbf_ao_mxnbf_ce*nbf_ao_mxnbf_ce
1230
if (.not.MA_Push_Get(MT_Dbl,nscr,'scr',lscr,iscr))
1231
& call errquit('dft_scf: cannot allocate scr',0, MA_ERR)
1232
call util_ga_mat_reduce(nbf_ao, natoms, int_mb(icetobfr),
1233
& g_dens, ipol, dbl_mb(irdens_atom),
1234
& 'rms', dbl_mb(iscr), nbf_ao_mxnbf_ce,.true.)
1235
if (.not.ma_pop_stack(lscr))
1236
& call errquit('dft_scf: cannot pop stacklscr:',0, MA_ERR)
1238
c Compute the XC potential and energy.
1241
call ga_zero(g_vxc(1))
1243
call ga_zero(g_vxc(2))
1246
if (cam_exch) call case_setflags(.true.) ! set attenuation
1248
if (oprint_time)call dft_tstamp(' Before call to GETVXC. ')
1250
& (rtdb, Exc, Ecoul,nExc, iVxc_opt, g_xcinv, g_dens,
1251
& g_vxc, IOLGC, g_wght, g_xyz, g_nq,wght_GA, rho_n,
1252
& dbl_mb(irdens_atom), int_mb(icetobfr), natoms)
1253
c write(*,*)"Ecoul, Exc(1)", Ecoul, Exc(1)
1256
if (cam_exch) call case_setflags(.false.) ! unset attenuation if set
1258
c == add in zora contributions ==
1261
c == calculate scalar energy ==
1262
Ezora_sf = ga_ddot(g_dens(1),g_zora_sf(1))
1263
& + ga_ddot(g_dens(2),g_zora_sf(2))
1264
Ecore = Ecore + Ezora_sf
1266
c == combine scalar part with the xc matrices ==
1267
call ga_dadd(1.d0,g_vxc(1),1.d0,g_zora_sf(1),g_vxc(1))
1268
call ga_dadd(1.d0,g_vxc(2),1.d0,g_zora_sf(2),g_vxc(2))
1270
c == add spin-orbit zora to fock matrix ==
1271
call ga_fock_so(g_zora_so, g_fockso, nbf_ao)
1275
c == add in the exchange-correlation to the fock matrix ==
1277
call ga_dadd_patch( 1.d0, g_fockso(1), 1, nbf_ao,
1279
& 1.0d0, g_vxc(1), 1, nbf_ao,
1281
& g_fockso(1), 1, nbf_ao,
1283
call ga_dadd_patch( 1.d0, g_fockso(1), 1+nbf_ao, nbf_mo,
1285
& 1.0d0, g_vxc(2), 1, nbf_ao,
1287
& g_fockso(1), 1+nbf_ao, nbf_mo,
1290
c == get the exact exchange contribution ==
1291
call xc_exso(rtdb,Exc,Ecoul,nExc,g_densso,g_fockso)
1294
& call dft_tstamp(' End of parallel region. ')
1296
c Calculate the total electronic energy.
1299
Etnew = Ecore + Ecoul + Exc(1)
1300
if(det_eng)goto 2001
1302
Etnew = Ecore + Ecoul + Exc(1) + Exc(2)
1303
if(det_eng)goto 2001
1306
if (last_time_energy)then
1308
c If open shell put the total density matrix back in
1309
c g_dens(1) and quit.
1311
call ga_dadd(one, g_dens(1), one, g_dens(2), g_dens(1))
1315
delta = abs(etold-etnew)
1320
homo_lumo_gap = 200.0d0
1322
c Symmetrize the Fock matrix
1325
& call sym_symmetrize(geom, AO_bas_han, .false., g_fock)
1327
call ga_symmetrize(g_fock)
1329
c DIIS step taken here.
1332
call diis_driver_so(toll_s, derr, icall, nfock,
1333
& nbf_mo, g_fockso, g_densso,
1334
& g_svecs, svals, diising, nodiis)
1338
g_scr = ga_create_atom_blocked(geom, AO_bas_han, 'ga scr')
1340
c Put s-1/2 in g_scr.
1343
call diis_bld12_so(toll_s, svals, g_svecs, g_scr,
1344
& g_tmp, nbf_ao, iw)
1346
c map s-1/2 to the nbf_mo by nbf_mo g_scr2
1348
if(.not.ga_create(mt_dbl, 2*nbf, 2*nbf,'scr2', 0, 0, g_scr2))
1349
& call errquit('dft_scf_so: error creating scr2',0, GA_ERR)
1350
call ga_zero(g_scr2)
1351
call ga_fock_sf(g_scr, g_scr2, nbf_ao)
1353
c Transform Fock matrix.
1355
call ga_zero(g_tmp_ri)
1356
call ga_dgemm('T', 'N', nbf_mo, nbf_mo, nbf_mo, one,
1357
& g_scr2, g_fockso(1), zero, g_tmp_ri)
1358
call ga_dgemm('N', 'N', nbf_mo, nbf_mo, nbf_mo, one,
1359
& g_tmp_ri, g_scr2, zero, g_fockso(1))
1361
call ga_zero(g_tmp_ri)
1362
call ga_dgemm('T', 'N', nbf_mo, nbf_mo, nbf_mo, one,
1363
& g_scr2, g_fockso(2), zero, g_tmp_ri)
1364
call ga_dgemm('N', 'N', nbf_mo, nbf_mo, nbf_mo, one,
1365
& g_tmp_ri, g_scr2, zero, g_fockso(2))
1367
c Level shifting is implemented here (similarity
1368
c transformation before standard eigensolver). Note,
1369
c levelshifting is appropriate once a transformation
1370
c is available which makes the resulting Fock matrix
1371
c diagonally dominant, e.g., in an approximate MO basis.
1372
c Also note, there are many matrix multiplies with S^+-1/2
1373
c which are redundant if one is sure that the former basis
1376
c levelshifting = .false.
1377
if (levelshifting)then
1379
c save the old vectors
1381
call ga_copy(g_moso(1), g_old(1))
1382
call ga_copy(g_moso(2), g_old(2))
1384
c Build a matrix which is diagonal in the "MO" rep,
1385
c back-transform, and shift the current Fock matrix
1387
c Use S^+1/2 * old movecs (as a transform).
1390
call diis_bld12_so(toll_s, svals, g_svecs, g_scr,
1391
& g_tmp, nbf_ao, iw)
1392
call ga_zero(g_scr2)
1393
call ga_fock_sf(g_scr, g_scr2, nbf_ao)
1394
call ga_zero(g_tmp_ri)
1395
call ga_dgemm('N', 'N', nbf_mo, nbf_mo, nbf_mo, one,
1396
& g_scr2, g_moso(1), zero, g_tmp_ri)
1397
call ga_copy(g_tmp_ri, g_moso(1))
1399
call ga_zero(g_tmp_ri)
1400
call ga_dgemm('N', 'N', nbf_mo, nbf_mo, nbf_mo, one,
1401
& g_scr2, g_moso(2), zero, g_tmp_ri)
1402
call ga_copy(g_tmp_ri, g_moso(2))
1404
c Build diagonal matrix.
1406
call ga_zero(g_tmp_ri)
1407
do j = nTotOcc+1+me, nbf_mo, nproc
1408
call ga_put(g_tmp_ri, j, j, j, j, rlshift, 1)
1411
c Transform this into "AO" basis and add to current
1414
call ga_dgemm('N', 'N', nbf_mo, nbf_mo, nbf_mo, one,
1415
& g_moso(1), g_tmp_ri, zero, g_scr2)
1416
call ga_dgemm('N', 'T', nbf_mo, nbf_mo, nbf_mo, one,
1417
& g_scr2, g_moso(1), one, g_fockso(1))
1418
call ga_dgemm('N', 'N', nbf_mo, nbf_mo, nbf_mo, one,
1419
& g_moso(2), g_tmp_ri, zero, g_scr2)
1420
call ga_dgemm('N', 'T', nbf_mo, nbf_mo, nbf_mo, one,
1421
& g_scr2, g_moso(2), one, g_fockso(1))
1423
call ga_dgemm('N', 'N', nbf_mo, nbf_mo, nbf_mo, one,
1424
& g_moso(1), g_tmp_ri, zero, g_scr2)
1425
call ga_dgemm('N', 'T', nbf_mo, nbf_mo, nbf_mo, mone,
1426
& g_scr2, g_moso(2), one, g_fockso(2))
1427
call ga_dgemm('N', 'N', nbf_mo, nbf_mo, nbf_mo, one,
1428
& g_moso(2), g_tmp_ri, zero, g_scr2)
1429
call ga_dgemm('N', 'T', nbf_mo, nbf_mo, nbf_mo, one,
1430
& g_scr2, g_moso(1), one, g_fockso(2))
1435
c Solve for the eigenvalues and eigenvectors of the Hamiltonian.
1437
call ga_symmetrize(g_fock)
1438
if (oprint_intermediate_fock)then
1440
cso#if defined(PARALLEL_DIAG)
1441
cso call ga_diag_std(g_fock, g_tmp, Dbl_MB(k_eval(ispin)))
1443
cso call ga_diag_std_seq(g_fock, g_tmp, Dbl_MB(k_eval(ispin)))
1445
cso call ga_diag_compl(g_fockso(1), g_fockso(2), g_moso(1),
1446
cso & g_moso(2), Dbl_MB(k_eval(1)))
1447
c write(*,*)"compare"
1450
DCpl_mb(ia+(nbf_mo)*(i-1)+(j-1))=dcmplx(0.0, 0.0)
1454
call ga_get(g_fockso(1), 1,i, i,i, dbl_mb(ibuff),1)
1456
c write(*,*)(dbl_mb(ibuff+ijk), ijk=0,nbf_mo-1)
1458
DCpl_mb(ia+(nbf_mo)*(i-1)+(j-1))=
1459
= dcmplx(dbl_mb(ibuff+j-1),0d0)
1461
call ga_get(g_fockso(2), 1,i, i,i, dbl_mb(ibuff),1)
1463
DCpl_mb(ia+(nbf_mo)*(i-1)+(j-1))=
1464
$ DCpl_mb(ia+(nbf_mo)*(i-1)+(j-1))
1465
$ +dcmplx(0d0,dbl_mb(ibuff+j-1))
1468
call ga_zero(g_moso(1))
1469
call ga_zero(g_moso(2))
1470
call zheev( 'V', 'U', nbf_mo, DCpl_mb(ia), nbf_mo,
1471
$ Dbl_mb(k_eval(1)),
1472
$ DCpl_mb(iwork), LLWORK, Dbl_mb(irwork), INFO )
1475
dbl_mb(ibuff+j-1)=0.0d0
1476
dbl_mb(ibuff+j-1)=dble(DCpl_mb(ia+nbf_mo*(i-1)+(j-1)))
1479
call ga_put(g_moso(1),1,nbf_mo,i1,i1,dbl_mb(ibuff),1)
1480
trace = ddot(nbf_mo,dbl_mb(ibuff),1,dbl_mb(ibuff),1)
1482
dbl_mb(ibuff+j-1)=0.0d0
1484
$ dimag(dcmplx(DCpl_mb(ia+nbf_mo*(i-1)+(j-1))))
1487
call ga_put(g_moso(2),1,nbf_mo,i1,i1,dbl_mb(ibuff),1)
1488
trace = ddot(nbf_mo,dbl_mb(ibuff),1,dbl_mb(ibuff),1)
1490
c write(*,'("before transform")')
1491
c write(*,*)(Dbl_mb(k_eval(1)+i),i=0,nbf_mo-1)
1492
c call ga_print(g_moso(1))
1493
c call ga_print(g_moso(2))
1495
c Check HOMO/LUMO gap.
1497
homo = Dbl_MB(k_eval(1)+nTotEl-1)
1498
lumo = Dbl_MB(k_eval(1)+nTotEl)
1500
c If levelshifting then tidy up.
1502
if (levelshifting)then
1504
c Put S^-1/2 back in g_scr2 (use g_fock as temp scr).
1507
call diis_bld12_so(toll_s, svals, g_svecs, g_scr,
1508
& g_fock, nbf_ao, iw)
1509
call ga_zero(g_scr2)
1510
call ga_fock_sf(g_scr, g_scr2, nbf_ao)
1513
c Back-transform eigenvectors with S^-1/2.
1515
call ga_dgemm('N', 'N', nbf_mo, nbf_mo, nbf_mo, one,
1516
& g_scr2, g_moso(1), zero, g_fockso(1))
1517
call ga_dgemm('N', 'N', nbf_mo, nbf_mo, nbf_mo, one,
1518
& g_scr2, g_moso(2), zero, g_fockso(2))
1519
if (.not. ga_destroy(g_scr)) call errquit
1520
& ('dft_scf: could not destroy g_scr', 0, GA_ERR)
1521
if (.not. ga_destroy(g_scr2)) call errquit
1522
& ('dft_scf_so: could not destroy g_scr2', 0, GA_ERR)
1524
c Keep orbital ordering according to principle
1525
c of maximum overlap with previous iteration.
1528
. call dft_mxovl(ao_bas_han, dbl_mb(k_eval(ispin)),
1529
& g_tmp, g_movecs(ispin), g_s,g_fock,
1530
. noc,oprint_conv_details,homo,lumo)
1532
call ga_zero(g_moso(1))
1533
call ga_zero(g_moso(2))
1534
call ga_copy(g_fockso(1), g_moso(1))
1535
call ga_copy(g_fockso(2), g_moso(2))
1537
c determine homo-lumo gap
1539
homo_lumo_gap = min(homo_lumo_gap, (lumo-homo-rlshift))
1540
if (me.eq.0 .and. oprint_conv_details)
1541
& write(LuOut,4224)homo,lumo,rlshift, homo_lumo_gap
1542
4224 format(10x,' HOMO = ',f6.2,' LUMO = ',f6.2,
1543
& ' RLSHIFT = ',f6.2,' HL_GAP = ',f6.2)
1547
c Save previous density for convergence check.
1549
call ga_copy(g_dens(1), g_movecs(1))
1550
call ga_copy(g_dens(2), g_movecs(2))
1552
c symmetry adapt vectors?
1555
call scf_sym_adapt_so(ao_bas_han, g_moso,
1556
& oprint_syma, 2*nbf_ao, name,
1561
c save the old density matrix for damping
1563
call ga_copy(g_densso(1), g_damp_so(1))
1564
call ga_copy(g_densso(2), g_damp_so(2))
1566
c Form a new density matrix.
1569
call ga_zero(g_densso(1))
1570
call ga_zero(g_densso(2))
1572
switch_sclMO_so=0 ! FA=09-26-11 set OFF scaleMO_so
1574
c the fractionally occupied mo's are scaled by the sqrt of the fon
1577
c scale = sqrt(avg_fon)
1578
do i = ncore_fon + 1, ncore_fon + nmo_fon
1579
if (i> nbf_mo) call errquit(
1580
& 'dft_densm:fon g_moso index exceeds nbf_mo',
1582
scale = sqrt(dbl_mb(i-1+k_occ)) ! sqrt(occupation number)
1583
call ga_scale_patch(g_moso(1),
1584
& 1, nbf_mo, i, i, scale)
1585
call ga_scale_patch(g_moso(2),
1586
& 1, nbf_mo, i, i, scale)
1588
if(me.eq.0) write(luout,'(5x,a)') 'FON applied'
1590
c ---- FA-02-15-11 : occupations keyword ---- START
1591
c call dft_scaleMO_so(rtdb,g_moso,dbl_mb(k_occ),g_densso,
1592
c & nbf_mo,nTotOcc,switch_sclMO_so)
1593
c ---- FA-02-15-11 : occupations keyword ---- END
1596
c if (switch_sclMO_so.ne.1) then
1597
c if (ga_nodeid().eq.0)
1598
c & write(*,*) 'ENTER calc so-density matrix std-way'
1599
call dft_densm_so(g_densso, g_moso, nbf_ao, nTotOcc)
1602
c restore the scaled mo's
1605
c scale = sqrt(avg_fon)
1606
do i = ncore_fon + 1, ncore_fon + nmo_fon
1607
if (i> nbf_mo) call errquit(
1608
& 'dft_densm:fon g_moso index exceeds nbf_mo',
1610
if (dbl_mb(i-1+k_occ) < 1d-4) call errquit(
1611
& 'dft_densm:fon frac occup < 1E-4. Aborting suspisciously',
1613
scale = 1d0/sqrt(dbl_mb(i-1+k_occ)) ! 1/sqrt(occupation number)
1614
call ga_scale_patch(g_moso(1),
1615
& 1, nbf_mo, i, i, scale)
1616
call ga_scale_patch(g_moso(2),
1617
& 1, nbf_mo, i, i, scale)
1621
cso g_dens(1)=Re(Daa)+Re(Dbb) and g_dens(2)=Re(Dbb)
1622
cso For coulomb and xc potentials only the alpha, Re(Daa), and
1623
cso beta, Re(Dbb) densities are needed
1625
call ga_zero(g_dens(1))
1626
call ga_zero(g_dens(2))
1627
call ga_dens_sf(g_dens, g_densso(1), nbf_ao)
1629
c check g_dens, calculate tr[P S] and print
1630
c if (debug_fon) call dft_pstrace(g_dens,ao_bas_han,nbf_ao,oskel)
1632
pstrace=ga_ddot(g_dens(1),g_s)
1633
pstrace=pstrace + ga_ddot(g_dens(2),g_s)
1634
if(ga_nodeid().eq.0) write (luout,'(5x,a,1x,e15.7)')
1635
& 'tr(P*S): ',pstrace
1636
if (.not. rtdb_put(rtdb, 'dft:pstrace', mt_dbl, 1, pstrace))
1637
& call errquit('dft_scf: rtdb_put pstrace failed', 1, RTDB_ERR)
1642
c Check convergence on Density.
1644
rms(1) = dft_dencvg(g_dens(1), g_movecs(1), nbf_ao)
1645
rms(2) = dft_dencvg(g_dens(2), g_movecs(2), nbf_ao)
1648
if (oprint_conv.and.iter.eq.1.and.me.eq.0)then
1649
nheap = MA_Inquire_Heap(MT_Dbl)
1650
nstack = MA_Inquire_Stack(MT_Dbl)
1652
write(LuOut,'(10x,a,f10.2,i20)')
1653
& ' Heap Space remaining (MW): ',dble(nheap)*1.D-06,nheap
1654
write(LuOut,'(10x,a,f10.2,i20)')
1655
& 'Stack Space remaining (MW): ',dble(nstack)*1.D-06,nstack
1656
call util_flush(LuOut)
1659
21 format(/,10x,' Memory utilization after 1st SCF pass: ')
1661
& 1x,' convergence iter energy DeltaE ',
1662
& 'RMS-Dens Diis-err time'/
1663
& 1x,'---------------- ----- ----------------- --------- ',
1664
& '--------- --------- ------')
1665
if (oprint_conv.and.me.eq.0)then
1666
current_cpu = util_cpusec()
1668
write(LuOut,2)ndamp,rlshift,
1670
& -etold+etnew,sqrt(rms(1)),derr(1),current_cpu
1671
if (ipol.eq.2)write(LuOut,3)sqrt(rms(2)),derr(2)
1673
write(LuOut,22)ndamp,rlshift,
1675
& -etold+etnew,sqrt(rms(1)), current_cpu
1676
if (ipol.eq.2)write(LuOut,23)sqrt(rms(2))
1678
call util_flush(LuOut)
1680
2 format(1x,'d=',i2,',ls=',f3.1,',diis',1x,i5,f18.10,
1681
& 1p,3d10.2,0p,f8.1)
1682
3 format(51x,1p,2d10.2)
1683
22 format(1x,'d=',i2,',ls=',f3.1,6x,i5,f18.10,
1684
& 1p,2d10.2,10x,0p,f8.1)
1685
23 format(51x,1p,1d10.2)
1689
call ecce_print1 ('iteration counter', mt_int, iter, 1)
1690
call ecce_print1 ('iterative total energy difference',
1691
& mt_dbl, -etold+etnew, 1)
1692
call ecce_print1 ('iterative total density difference',
1693
& mt_dbl, sqrt(rms(1)), 1)
1697
c save eigenvectors to movecs file
1699
if (.not.movecs_write_so
1700
$ (rtdb, ao_bas_han, movecs_out, 'sodft', title,
1701
& nbf_mo, dbl_mb(k_occ), dbl_mb(k_eval(1)), g_moso))
1702
& call errquit('dft_scf_so: movec_write failed', 0, DISK_ERR)
1706
if (me .eq. 0.and.oprint_eval)then
1707
if (util_print('intermediate evals', print_default))then
1708
call util_print_centered(LuOut,'eigenvalues',
1710
call output(dbl_mb(k_eval(1)), 1, min(nTotEl+10,nbf_mo),
1711
& 1, 1, nbf_mo, 1, 1)
1715
if (oprint_vecs)then
1718
call util_print_centered(LuOut,
1719
& 'Intermediate MO vectors',40,.true.)
1721
call util_flush(LuOut)
1725
c If open shell compute overlap of alpha orbitals with beta
1728
if ((ipol.gt.1).and.(oprint_interm_overlap)) then
1729
call dft_mxspin_ovlp(nbf_ao,nmo,ao_bas_han,g_movecs(1),
1730
& g_movecs(2), g_tmp)
1733
c computation of <S2> for open shell
1735
if ((ipol.gt.1).and.(oprint_interm_S2)) then
1737
call dft_s2_value(geom, AO_bas_han, .false., noc(1), noc(2),
1738
& nbf_ao, g_dens(1), g_dens(2))
1741
c Form the total density matrix.
1743
call ga_dadd(one, g_dens(1), one, g_dens(2), g_dens(1))
1746
c Check for SCF convergence.
1750
call dft_zora_scfcvg(rms, derr, Etold, Etnew,
1751
& e_conv, d_conv, g_conv, ipol,
1752
& iter, iterations, idone, rtdb,
1753
& converged, diising)
1755
call dft_scfcvg(rms, derr, Etold, Etnew,
1756
& e_conv, d_conv, g_conv, ipol,
1757
& iter, iterations, idone, rtdb,
1758
& converged, diising)
1760
if (delta.lt.1.d-3)then
1762
c Set coulomb acc to max (e.g., input parameter).
1763
c (note, may also require re-initializing DIIS)
1767
tol_rho = tol_rho_max
1771
c Damping implemented here.
1774
pp = dble(ndamp)*1.d-2
1776
call ga_dadd(pp, g_damp_so(1), onempp, g_densso(1),
1778
call ga_dadd(pp, g_damp_so(2), onempp, g_densso(2),
1780
call ga_zero(g_dens(1))
1781
call ga_zero(g_dens(2))
1782
call ga_dens_sf(g_dens, g_densso(1), nbf_ao)
1783
call ga_dadd(one, g_dens(1), one, g_dens(2), g_dens(1))
1791
c Check convergence parameters.
1793
if ((delta.lt.dampon.and.delta.gt.dampoff).or.iter.le.ncydp)then
1801
if ((delta.lt.levlon.and.delta.gt.levloff).or.
1802
& (iter.le.ncysh))then
1803
if (homo_lumo_gap.lt.hl_tol)then
1804
levelshifting = .true.
1805
rlshift = rlshift_input
1806
if (check_shift)then
1807
if (lumo .lt. homo)then
1808
levelshifting = .false.
1809
if (me.eq.0 .and. oprint_conv_details)
1810
& write(LuOut,2224)homo, lumo
1814
levelshifting = .false.
1815
rlshift = rlshift_def
1818
levelshifting = .false.
1819
rlshift = rlshift_def
1821
2224 format(10x,'HOMO = ',f6.2,' LUMO (with shift) = ',f6.2,
1822
& /,10x,'Unshifted LUMO is less than HOMO.',
1823
& /,10x,'Turning levelshifting OFF this iteration.')
1825
if ((delta.lt.diison.and.delta.gt.diisoff).or.
1826
& iter.le.ncyds.or.keep_diis_on)then
1829
c Once started, keep DIIS on until diisoff threshold.
1831
keep_diis_on = .true.
1835
if (delta.lt.diisoff.or.(ncyds.gt.0.and.iter.gt.ncyds))then
1837
keep_diis_on = .false.
1840
if (nodamping)damping = .false.
1841
if (nolevelshifting) then
1842
levelshifting = .false.
1845
if (nodiis)diising = .false.
1849
c activate disp if is present in rtdb
1850
c or if dispersion is automatically included in the xc functional
1851
if (.not.rtdb_get(rtdb, 'dft:disp', mt_log, 1, disp))
1854
if(disp.or.xc_chkdispauto())
1855
& call xc_vdw(rtdb,geom,Etnew,dum,'energy')
1860
if ((lumo - homo).lt.-hl_tol.and.lmaxov)then
1862
if (me.eq.0 .and. oprint_conv_details)
1863
& write(LuOut,224)homo, lumo
1864
224 format(10x,' HOMO = ',f6.2,' LUMO = ',f6.2,
1865
& /,10x,'Significant orbital reordering with',
1866
& ' maximum overlap',
1867
& /,10x,'turned ON. Turning max_ovl OFF.')
1870
if (oprint_energy_step.and.me.eq.0)then
1871
current_cpu = util_cpusec()
1873
write(LuOut,222)etnew+enuc, ecore, Ecoul, Exc(1), enuc,
1874
& rho_n, current_cpu
1876
write(LuOut,223)etnew+enuc, ecore, Ecoul, Exc(1), Exc(2),
1877
& enuc, rho_n, current_cpu
1881
c Check for remaining time to exit "gracefully"
1883
current_wall = util_wallsec()
1884
if ((iter-1).gt.1)then
1885
elapsed_wall = current_wall - save_wall
1886
save_wall = current_wall
1888
elapsed_wall = current_wall - start_wall
1889
save_wall = current_wall
1894
c If converged probably need a few seconds to clean things up
1895
c and calculate a few properties.
1897
wall_time_reqd = 5.0
1899
c == scale zora eigenvalues and energy ==
1902
call dft_zora_scale_so(
1903
& rtdb,g_dens_at,nexc, ! Added by FA
1914
& dbl_mb(k_eval(1)),
1924
c If not converged probably need at least the amount time
1925
c required for previous iteration (multiply by 1.2 to be on the safe side).
1927
wall_time_reqd = elapsed_wall*1.2d0
1929
int_wall_time_reqd = wall_time_reqd
1930
if (.not.util_test_time_remaining(rtdb, int_wall_time_reqd))then
1933
call util_print_centered(LuOut,
1934
& 'Exiting due to time limitations.', 20, .true.)
1939
if (idone.eq.0.or.(iswitc.lt.2.and.iter.lt.iterations))
1940
& go to 1000 ! begin new iteration
1941
if (idone.eq.1.and.(.not.last_time_energy))then
1942
last_time_energy = .true.
1943
go to 1000 ! build final total energies
1948
if (me.eq.0.and.oprint)then
1949
if (.not.converged)then
1951
call util_print_centered(LuOut,
1952
& 'Calculation failed to converge', 20, .true.)
1955
dft_time = dft_time+util_cpusec()
1958
write(LuOut,222)etnew+enuc,
1964
write(LuOut,223)etnew+enuc,
1971
if (do_zora) write(luout,2221) ener_scal
1972
write(luout,2222) rho_n
1973
write(luout,2223) dft_time
1976
& ' Total DFT energy =', f20.12/
1977
& ' One electron energy =', f20.12/
1978
& ' Coulomb energy =', f20.12/
1979
& ' Exchange-Corr. energy =', f20.12/
1980
& ' Nuclear repulsion energy =', f20.12/)
1983
& ' Total DFT energy =', f20.12/
1984
& ' One electron energy =', f20.12/
1985
& ' Coulomb energy =', f20.12/
1986
& ' Exchange energy =', f20.12/
1987
& ' Correlation energy =', f20.12/
1988
& ' Nuclear repulsion energy =', f20.12/)
1990
2221 format(' Scaling correction =', f20.12/)
1991
2222 format(' Numeric. integr. density =', f20.12/)
1992
2223 format(' Total iterative time =', f9.1,'s'//)
1994
call util_flush(LuOut)
1997
c print out the determinantal energies
2000
if (me.eq.0.and.oprint.and.det_eng)then
2002
c call util_print_centered(LuOut,
2003
c & 'Calculation failed to converge', 20, .true.)
2005
dft_time = dft_time+util_cpusec()
2007
write(LuOut,232)etnew+enuc, ecore, Ecoul, Exc(1), enuc,
2009
if (do_zora) write(luout,2221) ener_scal
2011
write(LuOut,233)etnew+enuc, ecore, Ecoul, Exc(1), Exc(2),
2012
& enuc, rho_n, dft_time
2013
if (do_zora) write(luout,2221) ener_scal
2016
& ' Determinant Energy'/
2017
& ' Total DFT energy =', f20.12/
2018
& ' One electron energy =', f20.12/
2019
& ' Coulomb energy =', f20.12/
2020
& ' Exchange-Corr. energy =', f20.12/
2021
& ' Nuclear repulsion energy =', f20.12//
2022
& ' Numeric. integr. density =', f20.12//
2023
& ' Total iterative time =', f9.1,'s'//)
2025
& ' Determinant Energy'/
2026
& ' Total DFT energy =', f20.12/
2027
& ' One electron energy =', f20.12/
2028
& ' Coulomb energy =', f20.12/
2029
& ' Exchange energy =', f20.12/
2030
& ' Correlation energy =', f20.12/
2031
& ' Nuclear repulsion energy =', f20.12//
2032
& ' Numeric. integr. density =', f20.12//
2033
& ' Total iterative time =', f9.1,'s'//)
2034
call util_flush(LuOut)
2037
c calculate the determinantal energy
2039
if(.not.fon) goto 2002
2040
if (.not.rtdb_get(rtdb,'sodft:ndet',mt_int,1,ndet)) ndet = 0
2042
if(idet .eq. 0 .and. ndet. gt. 0)then
2043
if(.not.ma_push_get(mt_int,nmo_fon*ndet,'det',kfon_occ,
2045
& call errquit('cannot alloctate lfon_occ',0, MA_ERR)
2046
if (.not.rtdb_get(rtdb,'sodft:occupancy',mt_int,
2047
& nmo_fon*ndet,int_mb(lfon_occ)))
2048
& call errquit('no occupancy specified',0, RTDB_ERR)
2049
endif ! (idet.eq.0) and (ndet.gt.0)
2051
if(idet.eq.ndet) goto 2002 ! all determinant energies have been calculated
2053
c Switch the mo order to match the determinantal occupancy
2055
iswap = ncore_fon + 1 ! first partially filled orbital
2056
jswap = iswap + idet ! forms the pair (i,j) to be swapped
2057
swap(1) = iswap ! assign the pair
2058
swap(2) = jswap ! assign the pair
2059
if (jswap.ge.iswap) then
2060
if(.not.rtdb_put(rtdb,'sodft:swap',mt_int,2,swap)) ! swap is queried in movecs_swap_so
2061
& call errquit('swap: failed to put nelem in rtdb', 0,
2063
call movecs_swap_so(rtdb,'dft',scftype,g_moso, ! do the swap
2064
& dbl_mb(k_occ),dbl_mb(k_eval(1)))
2067
c form a new density matrix and calculate the corresponding energy
2069
call dft_densm_so(g_densso,g_moso,nbf_ao,nTotEl)
2070
call ga_zero(g_dens(1))
2071
call ga_zero(g_dens(2))
2072
call ga_dens_sf(g_dens, g_densso(1), nbf_ao)
2073
call ga_dadd(one, g_dens(1), one, g_dens(2), g_dens(1))
2077
c Restore the mo order so that we have a fixed reference
2079
if (me.eq.0) write(luout,4443)
2080
4443 format(/,1x,'Restoring orbital order')
2082
if (jswap.ge.iswap) then
2083
if(.not.rtdb_put(rtdb,'sodft:swap',mt_int,2,swap)) ! swap is queried in movecs_swap_so
2084
& call errquit('swap: failed to put nelem in rtdb', 0,
2086
call movecs_swap_so(rtdb,'dft',scftype,g_moso, ! do the swap
2087
& dbl_mb(k_occ),dbl_mb(k_eval(1)))
2090
c Write determinant configuration
2094
& (int_mb(lfon_occ+idet*(nmo_fon)+i-1), i=1,nmo_fon)
2096
idet = idet + 1 ! determinant counter
2098
4444 format(/,1x,'Determinant Occupancy: ',8(1x,i3,1x))
2100
if(fon.and.(ndet.gt.0))then
2101
if (.not.ma_pop_stack(kfon_occ))
2102
& call errquit('dft_scf: cannot pop stack:lfon_occ',0,MA_ERR)
2107
if (.not. ga_destroy(g_damp_so(1))) call errquit
2108
& ('dft_scf_so: could not destroy g_damp_so', 0, GA_ERR)
2109
if (.not. ga_destroy(g_damp_so(2))) call errquit
2110
& ('dft_scf_so: could not destroy g_damp_so', 0, GA_ERR)
2112
c symmetry adapt vectors last time print symmetries, etc.
2115
c call scf_movecs_sym_adapt(ao_bas_han, g_movecs, oprint,
2116
c & nbf_ao, '- alpha', .true.,
2119
c & call scf_movecs_sym_adapt(ao_bas_han, g_movecs(2), oprint,
2120
c & nbf_ao, '- beta', .true.,
2121
c & int_mb(k_ir+nbf_ao))
2124
c Vector analysis stolen from rohf.F
2126
if (util_print('final vectors analysis', print_default)) then
2127
do ilo = 1,max(1,nTotEl-10)
2128
if (dbl_mb(k_eval(1)+ilo-1) .ge. eval_pr_tol_lo)
2131
961 do ihi = min(nTotEl+10,nbf_mo), nbf_mo
2132
if (dbl_mb(k_eval(1)+ihi-1) .ge. eval_pr_tol_hi)
2137
if (util_print('final vectors analysis', print_high)) then
2141
call movecs_anal_so(ao_bas_han, ilo, ihi, 0.15d0,
2143
& 'DFT Final Molecular Orbital Analysis',
2144
& .true., dbl_mb(k_eval(1)), oadapt,
2145
& int_mb(k_ir), .true., dbl_mb(k_occ))
2148
c call to Mulliken Pop Ananlysis
2153
& (' Total Density - Mulliken Population Analysis')
2154
call mull_pop(geom,ao_bas_han,g_dens(1),g_s, 'total')
2157
c analysis of spin density
2159
if (me.eq.0)call dft_header
2160
& (' Spin Density - Mulliken Population Analysis')
2161
call ga_dadd(one,g_dens(1),-2.d0,g_dens(2),g_dens(2))
2162
call mull_pop(geom,ao_bas_han,g_dens(2),g_s,'spin')
2163
call ga_dadd(one,g_dens(1),-1.d0,g_dens(2),g_dens(2))
2164
call ga_dscal(g_dens(2),0.5d0)
2168
c end infinite loop for SCF iterations
2170
c Store energy and convergence status ... must store before
2171
c write movecs since date of insertion is used.
2173
if (.not. rtdb_put(rtdb, 'sodft:energy',MT_DBL,1,(Etnew+Enuc)))
2174
& call errquit('dft_scf: failed to store energy in rtdb', 0,
2176
if (.not. rtdb_put(rtdb, 'sodft:converged',MT_LOG,1,converged))
2177
& call errquit('dft_scf: failed to store converged in rtdb',0,
2179
c if (rtdb_get(rtdb, 'sodft:converged', mt_log, 1, oconverged))
2180
c & write(*,*)"converged=", converged
2182
c output energies and eigenvectors to disk
2184
if (.not.movecs_write_so
2185
$ (rtdb, ao_bas_han, movecs_out, 'sodft', title,
2186
& nbf_mo, dbl_mb(k_occ), dbl_mb(k_eval(1)), g_moso))
2187
& call errquit('dft_scf_so: movec_write failed', 0, DISK_ERR)
2192
if (icall(1).gt.0)then
2194
call diis_driver_so(toll_s, derr, icall, nfock,
2195
& nbf_mo, g_fockso, g_densso,
2196
& g_svecs, svals, diising, nodiis)
2199
c If open shell compute overlap of alpha orbitals with beta orbitals.
2202
call dft_mxspin_ovlp(nbf_ao,nmo,ao_bas_han,g_movecs(1),
2203
& g_movecs(2), g_tmp)
2207
if (.not. ga_destroy(g_wght)) call errquit
2208
& ('dft_scf: could not destroy g_wght', 0, GA_ERR)
2209
if (.not. ga_destroy(g_xyz)) call errquit
2210
& ('dft_scf: could not destroy g_xyz', 0, GA_ERR)
2211
if (.not. ga_destroy(g_nq)) call errquit
2212
& ('dft_scf: could not destroy g_nq', 0, GA_ERR)
2215
c Restore alpha and beta densities.
2218
& call ga_dadd(one,g_dens(1),onem,g_dens(2),g_dens(1))
2220
c computation of <S2> for open shell
2223
call dft_s2_value(geom,AO_bas_han,.false.,noc(1),noc(2),
2224
& nbf_ao,g_dens(1),g_dens(2))
2227
c computation of multipole moments
2230
& call dft_mpole(rtdb, ao_bas_han, ipol, g_dens(1), g_dens(2))
2232
c print stolen for uhf.F
2234
if (util_print('schwarz',print_high).and.(.not.CDFIT))then
2235
call schwarz_print(natoms, nshells_ao)
2239
if (util_print('final evals', print_high))then
2240
call util_print_centered(LuOut,'Final eigenvalues',
2242
call util_print_centered(LuOut,
2243
& '(all occupied plus 10 virtual)',20,.true.)
2244
call output(dbl_mb(k_eval(1)),
2245
& 1, min(noc(1)+10,nbf_ao),
2246
& 1, 1, nbf_ao, 1, 1)
2250
if (oprint_final_vecs)then
2252
call util_print_centered(
2253
& LuOut,'Final MO vectors',40,.true.)
2255
call util_flush(LuOut)
2259
if (oprint_final_vecs)then
2260
cso call ga_print(g_movecs)
2261
cso if (ipol.eq.2)call ga_print(g_movecs(2))
2266
call movecs_ecce(nbf_ao, nmo, 1, nmo(1), dbl_mb(k_eval(1)),
2267
& dbl_mb(k_occ), int_mb(k_ir),
2268
& g_movecs(1), 'dft', 'alpha')
2269
if (ipol.eq.2)then ! spin-unrestricted
2270
call movecs_ecce(nbf_ao, nmo, 1, nmo(2), dbl_mb(k_eval(2)),
2271
& dbl_mb(k_occ+nbf_ao), int_mb(k_ir+nbf_ao),
2272
& g_movecs(2), 'dft', 'beta')
2274
call ecce_print1 ('total energy', mt_dbl, (Etold+Enuc), 1)
2275
call ecce_print1 ('nuclear repulsion energy', mt_dbl, Enuc, 1)
2276
call ecce_print1 ('coulomb energy', mt_dbl, Ecoul, 1)
2277
call ecce_print1 ('exchange energy', mt_dbl, Exc(1), 1)
2279
call ecce_print1 ('correlation energy', mt_dbl, Exc(2), 1)
2282
if (.not.ma_chop_stack(lbuff)) then
2283
call ma_summarize_allocated_blocks()
2286
call errquit('dft_scf: cannot pop stack:lbuff',lbuff, MA_ERR)
2289
if (.not.ma_pop_stack(lrwork))
2290
& call errquit('dft_scf: cannot pop stacklrwork',0, MA_ERR)
2291
if (.not.ma_pop_stack(lwork))
2292
& call errquit('dft_scf: cannot pop stack:lwork',0, MA_ERR)
2293
if (.not.ma_pop_stack(lw))
2294
& call errquit('dft_scf: cannot pop stack:lw',0, MA_ERR)
2295
if (.not.ma_pop_stack(la))
2297
& call errquit('dft_scf: cannot pop stack:la',0, MA_ERR)
2298
if (.not.ma_pop_stack(l_ir))
2299
& call errquit('dft_scf: cannot pop stack:l_ir',0, MA_ERR)
2302
if (.not. ga_destroy(g_fockt)) call errquit
2303
& ('dft_scf: could not destroy g_fockt', 0, GA_ERR)
2305
if (.not. ga_destroy(g_tmp)) call errquit
2306
& ('dft_scf: could not destroy g_tmp', 0, GA_ERR)
2308
call fock_2e_tidy(rtdb)
2311
call ecce_print_module_exit('dft', 'ok')
2313
call ecce_print_module_exit('dft', 'failed')
2316
c eval deallocation moved here from inside iteration loop
2318
if (.not.ma_pop_stack(l_eval))
2319
& call errquit('dft_scf: cannot pop stack:l_eval',0, MA_ERR)
2321
if (.not.ma_pop_stack(l_3cwhat))
2322
& call errquit('dft_scf: cannot pop stack:l_3cwhat',0, MA_ERR)
2323
if (.not.ma_pop_stack(l_3cERI))
2324
& call errquit('dft_scf: cannot pop stack:l_3cERI',0, MA_ERR)
2326
if (.not.ma_pop_stack(l_occ))
2327
& call errquit('dft_scf: cannot pop stack:l_occ',0, MA_ERR)
2328
if (.not.ma_pop_stack(lrdens_atom))
2329
& call errquit('dft_scf: cannot pop stack:lrdens_atom',0, MA_ERR)
2330
if (.not.ma_pop_stack(lcetobfr))
2331
& call errquit('dft_scf: cannot pop stack:lcetobfr',0, MA_ERR)
2332
if (.not.ma_pop_stack(lcntobfr))
2333
& call errquit('dft_scf: cannot pop stack:lcntobfr',0, MA_ERR)
2334
if (.not.ma_pop_stack(lcntoce))
2335
& call errquit('dft_scf: cannot pop stack:lcntoce',0, MA_ERR)
2337
if(.not.ga_destroy(g_moso(1)))
2338
& call errquit('dft_scf_so: error destroy Movecs Re',0, GA_ERR)
2339
if(.not.ga_destroy(g_moso(2)))
2340
& call errquit('dft_scf_so: error destroy Movecs Im',0, GA_ERR)
2341
if(.not.ga_destroy(g_fockso(1)))
2342
& call errquit('dft_scf_so: error destroy Fock Re',0, GA_ERR)
2343
if(.not.ga_destroy(g_fockso(2)))
2344
& call errquit('dft_scf_so: error destroy Fock Im',0, GA_ERR)
2345
if(.not.ga_destroy(g_densso(1)))
2346
& call errquit('dft_scf_so: error destroy DenMx Re',0, GA_ERR)
2347
if(.not.ga_destroy(g_densso(2)))
2348
& call errquit('dft_scf_so: error destroy DenMx Im',0, GA_ERR)
2349
if(.not.ga_destroy(g_tmp_ri))
2350
& call errquit('dft_scf_so: error destroy old re',0, GA_ERR)
2351
if(.not.ga_destroy(g_old(1)))
2352
& call errquit('dft_scf_so: error destroy old im',0, GA_ERR)
2353
if(.not.ga_destroy(g_old(2)))
2354
& call errquit('dft_scf_so: error destroy Tmp ReIm',0, GA_ERR)
2355
if(.not.ga_destroy(g_so(1)))
2356
& call errquit('dft_scf_so: error destroy so z',0, GA_ERR)
2357
if(.not.ga_destroy(g_so(2)))
2358
& call errquit('dft_scf_so: error destroy so y',0, GA_ERR)
2359
if(.not.ga_destroy(g_so(3)))
2360
& call errquit('dft_scf_so: error destroy so x',0, GA_ERR)
2362
c == deallocate zora arrays ==
2365
c == spin-free parts ==
2366
if (.not. ga_destroy(g_zora_sf(1))) call errquit(
2367
& 'dft_scf_so: ga_destroy failed zora sf 1',0, GA_ERR)
2368
if (.not. ga_destroy(g_zora_sf(2))) call errquit(
2369
& 'dft_scf_so: ga_destroy failed zora sf 2',0, GA_ERR)
2371
if (.not. ga_destroy(g_zora_scale_sf(1))) call errquit(
2372
& 'dft_scf_so: ga_destroy failed zora scale sf 1',0, GA_ERR)
2373
if (.not. ga_destroy(g_zora_scale_sf(2))) call errquit(
2374
& 'dft_scf_so: ga_destroy failed zora scale sf 2',0, GA_ERR)
2376
c == spin-orbit parts ==
2377
if(.not.ga_destroy(g_zora_so(1))) call errquit(
2378
& 'dft_scf_so: ga_destroy failed zora so z',0, GA_ERR)
2379
if(.not.ga_destroy(g_zora_so(2))) call errquit(
2380
& 'dft_scf_so: ga_destroy failed zora so y',0, GA_ERR)
2381
if(.not.ga_destroy(g_zora_so(3))) call errquit(
2382
& 'dft_scf_so: ga_destroy failed zora so x',0, GA_ERR)
2384
if(.not.ga_destroy(g_zora_scale_so(1))) call errquit(
2385
& 'dft_scf_so: ga_destroy failed zora so scale z',0, GA_ERR)
2386
if(.not.ga_destroy(g_zora_scale_so(2))) call errquit(
2387
& 'dft_scf_so: ga_destroy failed zora so scale y',0, GA_ERR)
2388
if(.not.ga_destroy(g_zora_scale_so(3))) call errquit(
2389
& 'dft_scf_so: ga_destroy failed zora so scale x',0, GA_ERR)
2393
dft_scf_so = converged
2396
if (.not. rtdb_get(rtdb, 'bgj:poliz', mt_log,
2397
& 1, do_poliz)) then
2401
c write(*,*)'*** dft_scf: calling cphf_poliz'
2402
if (.not. cphf_poliz(rtdb)) ! Never executed.
2403
$ call errquit(' cphf_poliz: failed from dft_scf !',0,
2409
1111 format(15x,'Core Energy: ',f20.10)
added
removed
.pc/09_backported_6.1.1_fixes.patch/src/nwdft/so_dft/dft_scf_so.F
