← Back to branch summary

~ubuntu-branches/ubuntu/trusty/nwchem/trusty-proposed

« back to all changes in this revision

Viewing changes to src/nwdft/scf_dft/dft_densm.F

  • Committer: Package Import Robot
  • Author(s): Michael Banck, Daniel Leidert, Andreas Tille, Michael Banck
  • Date: 2013-07-04 12:14:55 UTC
  • mfrom: (1.1.2)
  • Revision ID: package-import@ubuntu.com-20130704121455-5tvsx2qabor3nrui
Tags: 6.3-1
http://bugs.debian.org/696361
* New upstream release.
* Fixes anisotropic properties (Closes: #696361).
* New features include:
  + Multi-reference coupled cluster (MRCC) approaches
  + Hybrid DFT calculations with short-range HF 
  + New density-functionals including Minnesota (M08, M11) and HSE hybrid
    functionals
  + X-ray absorption spectroscopy (XAS) with TDDFT
  + Analytical gradients for the COSMO solvation model
  + Transition densities from TDDFT 
  + DFT+U and Electron-Transfer (ET) methods for plane wave calculations
  + Exploitation of space group symmetry in plane wave geometry optimizations
  + Local density of states (LDOS) collective variable added to Metadynamics
  + Various new XC functionals added for plane wave calculations, including
    hybrid and range-corrected ones
  + Electric field gradients with relativistic corrections 
  + Nudged Elastic Band optimization method
  + Updated basis sets and ECPs 

[ Daniel Leidert ]
* debian/watch: Fixed.

[ Andreas Tille ]
* debian/upstream: References

[ Michael Banck ]
* debian/upstream (Name): New field.
* debian/patches/02_makefile_flags.patch: Refreshed.
* debian/patches/06_statfs_kfreebsd.patch: Likewise.
* debian/patches/07_ga_target_force_linux.patch: Likewise.
* debian/patches/05_avoid_inline_assembler.patch: Removed, no longer needed.
* debian/patches/09_backported_6.1.1_fixes.patch: Likewise.
* debian/control (Build-Depends): Added gfortran-4.7 and gcc-4.7.
* debian/patches/10_force_gcc-4.7.patch: New patch, explicitly sets
  gfortran-4.7 and gcc-4.7, fixes test suite hang with gcc-4.8 (Closes:
  #701328, #713262).
* debian/testsuite: Added tests for COSMO analytical gradients and MRCC.
* debian/rules (MRCC_METHODS): New variable, required to enable MRCC methods.

Show diffs side-by-side

added added

removed removed

Lines of Context:
src/nwdft/scf_dft/dft_densm.F
 
1
c
1
2
      subroutine dft_densm(g_dens,g_vecs,nbf,nmo,ipol,
2
3
     ,     geom,AO_bas_han,noc,ntotel,
3
4
     ,  evals,focc,ssmear,tdots,
4
5
     ,     iter,search0,
5
6
     .     fon, nel_fon,nmo_fon,ncore_fon,
6
 
     .                     spinset,
7
 
     .                     rtdb) ! FA : for dft_scaleMO()
8
 
C$Id: dft_densm.F 21171 2011-10-07 18:21:53Z jochen $
9
 
C FA-10-04-11: Adding changes for using fractional occupations
 
7
     .                     spinset, rtdb)
 
8
c
10
9
      implicit none
 
10
c
11
11
#include "errquit.fh"
12
12
#include "bas.fh"
13
13
#include "geom.fh"
18
18
#include "dftps.fh"
19
19
#include "util.fh"
20
20
#include "rtdb.fh" 
21
 
 
 
21
c
22
22
      integer geom
23
23
      integer AO_bas_han
24
24
      integer ipol                  ! no. of spin states
33
33
      double precision ssmear       ! smearing sigma
34
34
      double precision tdots       ! total energy
35
35
      logical search0,fon
36
 
      double precision nel_fon(2)
37
 
      integer ncore_fon(2), nmo_fon(2)
 
36
      double precision nel_fon(4)
 
37
      integer ncore_fon(2), nmo_fon(4)
38
38
      integer ilo,ihi,jlo,jhi,numi,numj,icount
39
 
 
40
39
c
41
40
      integer i,isp,g_tmp
42
41
      integer me, nproc
48
47
      double precision util_erfc
49
48
      external util_erfc
50
49
      integer nocsave(2)
51
 
      logical dft_checkdg
52
 
      external dft_checkdg
 
50
c
 
51
      logical dft_fon
 
52
      external dft_fon
 
53
c
53
54
      double precision undovl
54
55
      parameter(undovl=-20d0*2.3025d0)
55
56
      parameter (zero=0.d0,toll=1.d-9,one=1.d0,
56
57
     ,   kbau=1.d0,eps=1.d-4)
57
 
      external dft_scaleMO ! FA-added-02-22-11
 
58
c
58
59
      integer rtdb         ! FA-02-22-11, for occupations
59
60
      integer switch_focc  ! FA-02-22-11, for occupations
60
61
      logical status       ! FA-02-22-11, for occupations
 
62
      double precision anoc(2),new_ntotel
61
63
      logical debug_fon
62
 
 
 
64
      logical no_avg_fon
 
65
      integer nocinit
 
66
      save nocinit,anoc
 
67
      data nocinit /1/
63
68
c
 
69
c     preliminaries
64
70
      sqrtpi=sqrt(acos(-1d0))
65
71
      me=ga_nodeid()
66
72
      nproc=ga_nnodes()
67
73
      rhfuhf=2d0
68
74
      if(ipol.eq.2) rhfuhf=1d0
69
 
      if(ssmear.lt.toll.or.iter.lt.-1) then
70
 
c       
71
 
c       check degener
72
 
c       
73
 
        
 
75
c       
 
76
c     check degener
 
77
c       
74
78
        if (rtdb_get(rtdb, 'dft:debugfon', mt_log, 1,
75
79
     &     debug_fon)) then
76
80
          if (debug_fon .and. me.eq.0)
77
 
     &       write(luout,*) "fon is on ",fon
 
81
     &       write(luout,*) "fractional occupation (fon) is on ",fon
78
82
        end if
79
 
                
 
83
c
 
84
c       average fractional occupations (leading orbitals)
80
85
        if (fon) then
81
 
          status = dft_checkdg(rtdb,nmo_fon,ncore_fon,nel_fon,
 
86
          status = dft_fon(rtdb,nmo_fon,ncore_fon,nel_fon,
82
87
     .       nbf,ntotel,focc,noc,ipol,me)
83
 
cng            if(dft_checkdg(nmo_fon,ncore_fon,nel_fon,
84
 
cng     .          nbf,ntotel,focc,noc,ipol,me)) then
85
88
          do isp=1,ipol
86
89
            call dft_focdm(focc(1+(isp-1)*nbf),noc(isp),geom,
87
 
     ,         AO_bas_han,nbf,
88
 
     *         g_vecs(isp),g_dens(isp),toll)
 
90
     &         AO_bas_han,nbf,g_vecs(isp),g_dens(isp),toll)
89
91
          enddo
90
 
          return
91
 
cng            endif
92
 
         endif
93
 
 
94
 
 
 
92
          if(ssmear.lt.toll)return
 
93
        endif ! fon
 
94
c
 
95
      if(ssmear.lt.toll.or.iter.lt.-1) then
 
96
c
95
97
        if (odftps) call pstat_on(ps_dgemm)
96
 
 
 
98
c
 
99
c       fractional occupation by orbital
97
100
        switch_focc=0
98
 
        status=rtdb_get(rtdb,'focc:occ-switch',
99
 
     &     mt_int,1,switch_focc)
 
101
        status=rtdb_get(rtdb,'focc:occup_switch', mt_int,1,switch_focc)
100
102
        if (switch_focc.eq.1 .and. status) then ! using specified occupations
101
103
          do isp=1,ipol
102
 
            call dft_scaleMO(rtdb,g_vecs, ! IN  : MO vectors
103
 
     &         focc,g_dens,     ! OUT : density matrix
104
 
     &         noc,             ! IN/OUT : orbital occupations
105
 
     &         nbf,ipol,ntotel)
106
 
          enddo                 ! end-loop-ispin
107
 
 
 
104
            call dft_frac_mo(rtdb,g_vecs(isp),focc,nbf,ipol,ntotel)
 
105
            call dft_focdm(focc(1+(isp-1)*nbf),noc(isp),geom,
 
106
     &         AO_bas_han,nbf,g_vecs(isp),g_dens(isp),toll)
 
107
          enddo  ! isp
 
108
          return
108
109
        else                    ! default occupations
109
110
          do isp=1,ipol
110
111
            call ga_dgemm('N', 'T', nbf, nbf, noc(isp),
111
112
     $         2d0/dble(ipol), g_vecs(isp),g_vecs(isp),
112
113
     $         zero, g_dens(isp))
113
114
          enddo
114
 
 
115
115
        endif                   ! switch_focc
116
116
 
117
117
        if (odftps) call pstat_off(ps_dgemm)
118
 
      else      
 
118
c
 
119
      else        ! with smearing
 
120
c
119
121
        if (.not. MA_Push_Get(MT_Dbl, nbf, 'tmpm', ltmpm, itmpm))
120
122
     &     call errquit('dftdensm: failed to alloc tmpm',nbf, MA_ERR)
 
123
c
 
124
c     fon
 
125
c
 
126
        if(nocinit.eq.1) then
 
127
           if(fon) then
 
128
              anoc(1)=ncore_fon(1)+nel_fon(1)
 
129
              anoc(2)=ncore_fon(2)+nel_fon(2)
 
130
           else
 
131
              anoc(1)=noc(1)
 
132
              anoc(2)=noc(2)
 
133
           endif
 
134
           nocinit=0
 
135
        endif
 
136
        new_ntotel=rhfuhf*(anoc(1)+anoc(2))
121
137
c       
122
138
c       initialize ef
123
139
c       
130
146
          call dfill(nbf*ipol, 0.d0, focc, 1)
131
147
          if(spinset.and.ipol.eq.2) then
132
148
            nmo(2)=0
133
 
            call dft_zero(2,nbf,nmo,noc(1),efermi(1),evals,
 
149
            call dft_zero(2,nbf,nmo,anoc(1),efermi(1),evals,
134
150
     ,         ssmear,toll,.true.)
135
 
            call dft_zero(2,nbf,nmo,noc(2),efermi(2),evals(nbf+1),
 
151
            call dft_zero(2,nbf,nmo,anoc(2),efermi(2),evals(nbf+1),
136
152
     .         ssmear,toll,.true.)
137
153
            nmo(2)=nmo(1)
138
154
          else
139
 
            call dft_zero(ipol,nbf,nmo,ntotel,efermi(1),evals,ssmear,
 
155
cold            call dft_zero(ipol,nbf,nmo,ntotel,efermi(1),evals,ssmear,
 
156
            call dft_zero(ipol,nbf,nmo,new_ntotel,
 
157
     E            efermi(1),evals,ssmear,
140
158
     .         toll,spinset)
141
159
            efermi(2)=efermi(1)
142
160
          endif
199
217
c       
200
218
c       compute entropy correction to total e
201
219
c     
202
 
        if(.not.fon) then
 
220
        if(ssmear.gt.1d-8) then
203
221
          cksum=0.d0
204
222
         do i=1,nmo(1)
205
223
           x=(evals(i)-efermi(1))/ssmear
223
241
C          1723       format(' Foccs ',(
224
242
         endif
225
243
       endif
 
244
c
226
245
       if(spinset.and.ipol.eq.2) then
227
246
         if(search0) then
228
247
           noc(1)=nint(ntot(1))
232
251
           noc(2)=nocsave(2)
233
252
         endif
234
253
       endif
235
 
      endif
 
254
c
 
255
      endif  ! smear
 
256
c
236
257
      return
237
258
      end
 
259
c
238
260
      subroutine dft_zero(ipol,nbf,nmo,
239
261
     ,     ntotel,efermi,evals,ssmear,toll,spinset)
240
262
      implicit none
247
269
      double precision evals(nbf*ipol) 
248
270
      double precision ssmear
249
271
      double precision toll
250
 
      integer ntotel
 
272
cold      integer ntotel
 
273
      double precision ntotel
251
274
      double precision efermi
252
275
      logical spinset
253
276
c
265
288
c
266
289
c     closed shell
267
290
c
268
 
         efermi=evals(ntotel/2)
 
291
         efermi=evals(int(ntotel/2))
269
292
      else
270
293
c
271
294
c     open shell
272
295
c
273
296
         if(spinset) then
274
 
            efermi=evals(ntotel/2)
 
297
            efermi=evals(int(ntotel/2))
275
298
         else
276
 
            efermi=0.5d0*(evals(ntotel/2)+evals(nbf+ntotel/2))
 
299
            efermi=0.5d0*(evals(int(ntotel/2))+
 
300
     +           evals(nbf+int(ntotel/2)))
277
301
         endif
278
302
      endif
279
303
      step=max(2.d0*ssmear,1.d0)
353
377
 
354
378
      return
355
379
      end
 
380
c
356
381
      double precision function fzero(
357
382
     *     ipol,efermi,evals,nbf,nmo,ntotel,ssmear)
358
383
      implicit none
360
385
      integer nbf,nmo(2)
361
386
      double precision efermi
362
387
      double precision evals(nbf*ipol)
363
 
      integer ntotel 
 
388
cold      integer ntotel 
 
389
      double precision ntotel
364
390
      double precision ssmear
365
391
#include "msgids.fh"
366
392
c
369
395
      external util_erfc
370
396
      double precision one,x,xb,rhfuhf,
371
397
     , fzeroa ,fzerob ,na,nb
372
 
      parameter(one=1.d0)
373
398
      rhfuhf=2d0/ipol
374
399
      fzeroa=0d0
375
400
      fzerob=0d0
384
409
         fzeroa=fzeroa+na
385
410
         if(ipol.eq.2) fzerob=fzerob+nb
386
411
      enddo
387
 
C      write(0,*) ' ef a b ',efermi,fzeroa,fzerob
388
 
C      write(0,*) 'abs  a b ',abs(noc(1)-fzeroa),abs(noc(2)-fzerob)
 
412
c      write(6,*) ' ef a b ',efermi,fzeroa,fzerob
389
413
         fzero=ntotel-(fzeroa+fzerob)*rhfuhf
390
 
C      write(0,*) ' fzero efer ',fzero,efermi
391
 
C      write(0,*) ' ==============='
 
414
c      write(6,*) ' ntotel  ',ntotel
 
415
c      write(6,*) ' fzero efer ',fzero,efermi
 
416
c      write(6,*) ' ==============='
392
417
      return
393
418
      end
 
419
c
394
420
      subroutine dft_focdm(focc,noc,geom,AO_bas_han,nbf,
395
421
     *     g_vecs,g_dens,toll)
396
422
      implicit none
408
434
      double precision zero
409
435
      parameter(zero=0.d0)
410
436
      logical debug
411
 
 
 
437
c
 
438
c     preliminaries
412
439
      me = ga_nodeid()
413
440
      nproc = ga_nnodes()
414
441
      nocold=noc
415
442
      debug = .false.
416
 
 
417
 
      if (me.eq.0 .and .debug) then
418
 
        write (luout,*) 'hello from dft_focdm'
419
 
      end if
420
 
        
 
443
c
421
444
      if (.not. MA_Push_Get(MT_Dbl, nbf, 'tmpm', ltmpm, itmpm))
422
445
     &   call errquit('dftdensm: failed to alloc tmpm',0, MA_ERR)
423
446
      if (.not. ga_duplicate(g_dens, g_tmp, 'ga_temp'))
444
467
c     reset noc for future calls
445
468
c     
446
469
      noc=nocold
447
 
c     call dfill(nbf, 0.d0, focc, 1)
448
 
c     call dfill(noc, rhfuhf, focc, 1)
449
470
      return
450
471
      end
451
 
 
452
 
      logical function dft_checkdg(rtdb,nmo_fon,ncore_fon,nel_fon,
453
 
     .     nbf,ntotel,focc,noc,ipol,
454
 
     , me)
455
 
c
456
 
c ... jochen 10/11: changed this routine to apply fractional occupations
457
 
c     during DFT SCF cycles. This also allows for the use of fractional
458
 
c     electron numbers overall.
459
 
c     This routine is called from within dft_densm during non-spinorbit
460
 
c     calculations
 
472
c
 
473
c     determine fractional occupations for fon calculations
 
474
c
 
475
      logical function dft_fon(rtdb,nmo_fon,ncore_fon,nel_fon,
 
476
     .     nbf,ntotel,focc,noc,ipol,me)
461
477
c
462
478
      implicit none
463
479
#include "global.fh"
465
481
#include "errquit.fh"
466
482
#include "rtdb.fh"
467
483
#include "mafdecls.fh"
 
484
c
468
485
      integer rtdb
469
486
      integer ipol
470
487
      integer nbf
474
491
      integer me
475
492
c
476
493
      integer i
477
 
      double precision avg_fon
478
 
      double precision nel_fon(2)
479
 
      integer nmo_fon(2), ncore_fon(2)
 
494
      double precision avg_fon, avg_fon2
 
495
      double precision nel_fon(4)
 
496
      integer nmo_fon(4), ncore_fon(2)
480
497
      integer ispin
481
498
      double precision ncheck
482
499
      logical debug
483
 
      integer ncor
484
 
 
 
500
      integer nfilled
 
501
      logical do_avg_fon
 
502
      logical do_core_fon
 
503
c
485
504
      debug = .false.
486
 
      if (rtdb_get(rtdb, 'dft:debugfon', mt_log, 1,
487
 
     &   debug)) continue
488
 
 
 
505
      if (.not. rtdb_get(rtdb, 'dft:debugfon', mt_log, 1, debug)) 
 
506
     &     debug = .false.
 
507
c
 
508
c     do average fractional occupation by default
 
509
      do_avg_fon = .true.
 
510
      if (.not.rtdb_get(rtdb,'dft:avg_fon',mt_log,1,do_avg_fon))
 
511
     &     do_avg_fon = .true.
 
512
c
 
513
c     do average occupation starting with core orbitals
 
514
      do_core_fon = .false.
 
515
      if (.not.rtdb_get(rtdb,'dft:core_fon',mt_log,1,do_core_fon))
 
516
     &     do_core_fon = .false.
 
517
c
489
518
      if (me.eq.0 .and. debug) then
490
519
        write (luout,*) 'FON: ipol, noc, ntotel',ipol,noc(:),ntotel
 
520
        write (luout,*) 'do_avg_fon: ', do_avg_fon
491
521
      end if                    ! debug
492
522
c
493
 
c  assignments      
494
 
c
495
 
c ... jochen: this functionality was not doing what 
496
 
c     I thought it was doing. Let's try differently
497
 
 
498
 
c$$$      do ispin=1,ipol
499
 
c$$$         avg_fon = dble(nel_fon(ispin))/dble(nmo_fon(ispin)) 
500
 
c$$$         do i=1,noc(ispin)-nel_fon(ispin)
501
 
c$$$            focc(i,ispin) = 2d0/ipol
502
 
c$$$         enddo
503
 
c$$$            ncheck=(noc(ispin)-nel_fon(ispin)+1d0)*
504
 
c$$$     .      2d0/dble(ipol)
505
 
c$$$         do i = noc(ispin)-nel_fon(ispin)+1,
506
 
c$$$     ,        noc(ispin)-nel_fon(ispin)+nmo_fon(ispin)
507
 
c$$$            focc(i,ispin) = avg_fon*(2d0/ipol)
508
 
c$$$            ncheck=ncheck+focc(i,ispin)
509
 
c$$$         enddo
510
 
c$$$      enddo
511
 
 
512
 
c ... jochen: new code:
513
 
 
514
523
      ncheck = 0d0
515
524
      do ispin = 1,ipol
516
 
 
517
 
        if (nmo_fon(ispin).lt.1) call errquit(
518
 
     &     'dft_densm:fon nmo_fon <1',
519
 
     &     1, INPUT_ERR)
520
 
        if (nel_fon(ispin).lt.0d0) call errquit(
521
 
     &     'dft_scf_so:fon nel_fon <0',
522
 
     &     1, INPUT_ERR)
523
 
 
524
 
        avg_fon = nel_fon(ispin)/dble(nmo_fon(ispin))
525
 
        ncor = ncore_fon(ispin)
526
 
        do i = 1,ncor
527
 
          if (i>nbf) call errquit(
528
 
     &       'dft_densm:fon focc index exceeds nbf',
529
 
     &       1, INPUT_ERR)
530
 
          focc(i,ispin) = 2d0/ipol
531
 
          ncheck = ncheck + focc(i,ispin)
532
 
        end do
533
 
        do i = 1,nmo_fon(ispin)   
534
 
          if (i+ncor>nbf) call errquit(
535
 
     &       'dft_densm:fon focc index exceeds nbf',
536
 
     &       1, INPUT_ERR)         
537
 
          focc(i+ncor,ispin) = avg_fon
538
 
          ncheck = ncheck + focc(i+ncor,ispin)
539
 
        end do
 
525
        if (nmo_fon(ispin).lt.1) 
 
526
     &     call errquit('dft_densm:fon nmo_fon <1',1, INPUT_ERR)
 
527
        if (nel_fon(ispin).lt.0d0) 
 
528
     &     call errquit('dft_scf_so:fon nel_fon <0',1, INPUT_ERR)
 
529
c
 
530
        avg_fon = nel_fon(ispin)/dble(nmo_fon(ispin))  ! average occupation
 
531
        if (nmo_fon(ispin+2).lt.1.or.nel_fon(ispin+2).le.0d0) then
 
532
            avg_fon2 = 0.0d0
 
533
        else
 
534
            avg_fon2 = nel_fon(ispin+2)/dble(nmo_fon(ispin+2))  ! average occupation
 
535
        endif
 
536
        nfilled = ncore_fon(ispin)  ! number of filled orbitals for each spin
 
537
c
 
538
c       initialize
 
539
        do i = 1,nbf
 
540
          focc(i,ispin) = 0.0
 
541
        end do   
 
542
c
 
543
c       fill the molecular orbitals either starting from the core or valence states
 
544
        if (do_core_fon) then
 
545
c
 
546
c         partially filled molecular orbitals
 
547
          do i = 1,nmo_fon(ispin)   
 
548
            if (i > nbf) call errquit(
 
549
     &       'dft_densm:fon focc index exceeds nbf',1,INPUT_ERR)
 
550
            if (do_avg_fon) then
 
551
             focc(i,ispin) = avg_fon  ! assign average occupation
 
552
            else
 
553
             focc(i,ispin) = nel_fon(ispin) ! assign given fractional electron
 
554
            end if ! do_avg_fon
 
555
            ncheck = ncheck + focc(i,ispin)
 
556
            if (me.eq.0 .and. debug) 
 
557
     &        write(luout,*) i,ispin,focc(i,ispin)
 
558
          end do  ! nmo_fon
 
559
c
 
560
c         fully filled molecular orbitals        
 
561
          do i = 1,nfilled
 
562
            if (i+nmo_fon(ispin)>nbf) call errquit(
 
563
     &       'dft_densm:fon focc index exceeds nbf',1,INPUT_ERR)
 
564
            focc(i+nmo_fon(ispin),ispin) = 2d0/ipol
 
565
            ncheck = ncheck + focc(i+nmo_fon(ispin),ispin)
 
566
            if (me.eq.0 .and. debug) 
 
567
     &        write(luout,*) i,ispin,focc(i+nmo_fon(ispin),ispin)
 
568
          end do  ! nfilled
 
569
 
 
570
        else
 
571
c
 
572
c         fully filled molecular orbitals
 
573
          do i = 1,nfilled
 
574
            if (i>nbf) call errquit(
 
575
     &       'dft_densm:fon focc index exceeds nbf',1,INPUT_ERR)
 
576
            focc(i,ispin) = 2d0/ipol
 
577
            ncheck = ncheck + focc(i,ispin)
 
578
            if (me.eq.0 .and. debug) 
 
579
     &        write(luout,*) i,ispin,focc(i,ispin)
 
580
          end do  ! nfilled
 
581
c
 
582
c         partially filled molecular orbitals
 
583
          do i = 1,nmo_fon(ispin)   
 
584
            if (i+nfilled > nbf) call errquit(
 
585
     &       'dft_densm:fon focc index exceeds nbf',1,INPUT_ERR)
 
586
            if (do_avg_fon) then
 
587
             focc(i+nfilled,ispin) = avg_fon  ! assign average occupation
 
588
            else
 
589
             focc(i+nfilled,ispin) = nel_fon(ispin) ! assign given fractional electron
 
590
            end if ! do_avg_fon
 
591
            ncheck = ncheck + focc(i+nfilled,ispin)
 
592
            if (me.eq.0 .and. debug) 
 
593
     &        write(luout,*) i,ispin,focc(i+nfilled,ispin)
 
594
          end do  ! nmo_fon
 
595
c
 
596
c         second partially filled molecular orbitals
 
597
          do i = 1,nmo_fon(ispin+2)
 
598
            if (i+nfilled > nbf) call errquit(
 
599
     &       'dft_densm:fon focc index exceeds nbf',1,INPUT_ERR)
 
600
            if (do_avg_fon) then
 
601
             focc(i+nfilled+nmo_fon(ispin),ispin) = avg_fon2  ! assign average occupation
 
602
            else
 
603
             focc(i+nfilled+nmo_fon(ispin),ispin) = nel_fon(ispin+2) ! assign given fractional electron
 
604
            end if ! do_avg_fon
 
605
            ncheck = ncheck + focc(i+nfilled+nmo_fon(ispin),ispin)
 
606
            if (me.eq.0 .and. debug)
 
607
     &       write(luout,*) i,ispin,focc(i+nfilled+nmo_fon(ispin),ispin)
 
608
          end do  ! nmo_fon second partial
 
609
       end if  ! do_core_fon
540
610
      end do                    ! ispin
541
 
      
 
611
c
542
612
      if (me.eq.0 .and. debug) then
543
 
        write (luout,*) 'FON: focc:',focc(:,1)
 
613
        write(luout,*) 'FON: focc:',focc(:,1)
 
614
        write(luout,*) 'ncheck:' , ncheck
 
615
        write(luout,*) 'ntotel:' , ntotel
544
616
      end if   
545
 
 
 
617
c
546
618
      if(abs(ncheck-dble(ntotel)).gt.1d-3 .and. me.eq.0) then
547
 
c         write(luout,*) ' frac. electrons ',ncheck,' vs ',ntotel
548
 
         dft_checkdg=.false.
 
619
         dft_fon=.false.
549
620
      else
550
 
         dft_checkdg=.true.
 
621
         dft_fon=.true.
551
622
      endif
552
 
c      if(me.eq.0) write(luout,'(5x,a)') 'FON applied'
 
623
c
553
624
      return
554
625
      end
 
626
c $Id: dft_densm.F 24151 2013-05-02 01:03:46Z edo $