← Back to branch summary

~siesta-pseudos-bases/siesta/trunk-psml

« back to all changes in this revision

Viewing changes to Src/m_new_dm.F90

  • Committer: Alberto Garcia
  • Date: 2019-09-02 14:09:43 UTC
  • mfrom: (427.6.323 trunk)
  • Revision ID: albertog@icmab.es-20190902140943-mzmbe1jacgefpgxw
Sync to trunk-776 (notably nc/soc wavefunction support)


Show diffs side-by-side

added added

removed removed

Lines of Context:
Src/m_new_dm.F90
67
67
    use m_ts_global_vars,only: TSrun
68
68
    use m_ts_electype, only : copy_DM
69
69
    use m_ts_options, only : TS_analyze
70
 
    use m_ts_options, only : N_Elec, Elecs, DM_bulk
 
70
    use m_ts_options, only : N_Elec, Elecs
71
71
    use m_ts_method
72
72
    use m_energies, only: Ef
73
73
 
79
79
    logical :: DM_init ! Initialize density matrix from file or from atomic density
80
80
    logical :: read_DM
81
81
    integer :: DM_in_history, n_depth
 
82
 
 
83
    ! Variable to signal how the initialization was performed.
 
84
    ! The following table lists the (current) possibilities:
 
85
    !  0: DM is filled from atomic information (no information present)
 
86
    !  1: .DM file is read, Siesta continuation
 
87
    !  2: An extrapolation of previous geometries.
 
88
    !     The DM's from the previous coordinates are kept in memory
 
89
    !  3: .TSDE file is read, TranSiesta continuation
 
90
    ! If the value is negative it means that one cannot mix in the
 
91
    ! first step, unless SCF.Mix.First.Force is true!
82
92
    integer :: init_method
83
93
 
84
94
    real(dp), pointer :: DM(:,:), EDM(:,:)
192
202
          call print_type(DM_2D)
193
203
       end if
194
204
 
195
 
       ! Defines the initialiazation to be re-using
196
 
       init_method = -1
 
205
       ! A negative value specifies that one cannot mix in the initial step
 
206
       init_method = -3
197
207
 
198
208
    end if
199
209
 
218
228
       ! If the Fermi-level has not been
219
229
       ! set, we initialize it to the mean of the
220
230
       ! electrode chemical potentials
221
 
       if ( DM_bulk > 0 ) then
 
231
       if ( any(Elecs(:)%DM_init > 0) ) then
222
232
 
223
233
          set_Ef = abs(Ef) < 0.00001_dp .and. &
224
234
               (.not. fdf_defined('TS.Fermi.Initial') ) 
244
254
            ! Not in buffer
245
255
            na_a = 0
246
256
            do iElec = 1 , na_u
247
 
              if ( a_isElec(iElec) ) na_a = na_a + 1
 
257
              if ( a_isElec(iElec) ) then
 
258
                if ( Elecs(atom_type(iElec))%DM_init > 0 ) then
 
259
                  na_a = na_a + 1
 
260
                end if
 
261
              end if
248
262
            end do
249
263
            allocate(allowed_a(na_a))
250
264
            na_a = 0 
251
265
            do iElec = 1 , na_u
252
266
              if ( a_isElec(iElec) ) then
253
 
                na_a = na_a + 1
254
 
                allowed_a(na_a) = iElec
 
267
                if ( Elecs(atom_type(iElec))%DM_init > 0 ) then
 
268
                  na_a = na_a + 1
 
269
                  allowed_a(na_a) = iElec
 
270
                end if
255
271
              end if
256
272
            end do
257
273
          end if
258
274
 
259
275
          do iElec = 1 , N_Elec
260
276
 
 
277
            ! We shift the mean by one fraction of the electrode
 
278
            if ( set_Ef ) then
 
279
              Ef = Ef + Elecs(iElec)%Ef / N_Elec
 
280
            end if
 
281
 
 
282
            if ( Elecs(iElec)%DM_init == 0 ) cycle
 
283
            
261
284
            if ( IONode ) then
262
285
              write(*,'(/,2a)') 'transiesta: Reading in electrode TSDE for ', &
263
286
                  trim(Elecs(iElec)%Name)
269
292
            call copy_DM(Elecs(iElec),na_u,xa,lasto,nsc,isc_off, &
270
293
                ucell, DM_2D, EDM_2D, na_a, allowed_a)
271
294
            
272
 
            ! We shift the mean by one fraction of the electrode
273
 
            if ( set_Ef ) then
274
 
              Ef = Ef + Elecs(iElec)%Ef / N_Elec
275
 
            end if
276
 
            
277
295
          end do
278
296
 
279
297
          ! Clean-up
290
308
 
291
309
          old_Ef = Ef
292
310
          Ef = fdf_get('TS.Fermi.Initial',Ef,'Ry')
293
 
          if ( init_method == 2 ) then
 
311
          if ( abs(init_method) == 3 ) then
294
312
             ! As the fermi-level has been read in from a previous
295
313
             ! calculation (TSDE), the EDM should only be shifted by the difference
296
314
             diff_Ef = Ef - old_Ef
301
319
 
302
320
          if ( IONode ) then
303
321
             write(*,*) ! new-line
304
 
             if ( init_method < 2 ) then
 
322
             if ( abs(init_method) < 3 ) then
305
323
                write(*,'(a,f9.5,a)')'transiesta: Setting the Fermi-level to: ', &
306
324
                     Ef / eV,' eV'
307
 
             else if ( init_method == 2 ) then
 
325
             else if ( abs(init_method) == 3 ) then
308
326
                write(*,'(a,2(f10.6,a))')'transiesta: Changing Fermi-level from -> to: ', &
309
327
                     old_Ef / eV,' -> ',Ef / eV, ' eV'
310
328
             end if
325
343
    ! Determine how the mixing of the first step should be performed
326
344
    !
327
345
    ! The idea is that one will allow mixing on the first SCF step
328
 
    ! only for atomicly filled orbitals. This will in most cases
 
346
    ! only for atomically filled orbitals. This will in most cases
329
347
    ! be a good initial choice, while in certain systems (spin-orbit)
330
348
    ! it may not be so good.
331
349
    ! Subsequently for any MD steps beyond the initial step, we will not
336
354
    ! iterations. For instance when performing FC runs (externally driven)
337
355
    ! this would lead to non-degenerate transverse eigenmodes for simple molecules
338
356
    if ( init_method == 0 ) then ! atomicly filled data
339
 
       ! allow mix_scf_first
340
 
    else if ( mix_scf_first ) then
 
357
      ! allow mix_scf_first
 
358
    else if ( mix_scf_first_force ) then
 
359
      ! user requested a mix in the first step
 
360
      mix_scf_first = .true.
 
361
      if ( IONode .and. init_method < 0 ) then
 
362
        ! Warn the user about this
 
363
        write(*,"(a)") "Mixing first scf iteration is *forced* although the sparsity pattern has changed..."
 
364
        write(*,"(a)") "Mixing non-compatible matrices might result in problems."
 
365
        write(*,"(a)") "Please use a linear mixer for the first few steps to remove history effects."
 
366
      end if
 
367
    else if ( mix_scf_first .and. init_method < 0 ) then
 
368
       ! Do not allow mixing first SCF step since we are reusing a DM
 
369
       ! from another geometry.
341
370
       mix_scf_first = .false.
342
371
    end if
343
372
    
385
414
    ! integer init_method           : returns method it has read the data with
386
415
    !                                   0 == atomic filling (possibly user-defined
387
416
    !                                   1 == .DM read
388
 
    !                                   2 == .TSDE read
 
417
    !                                   3 == .TSDE read
389
418
    ! *******************************************************************
390
419
 
391
420
    ! The spin-configuration that is used to determine the spin-order.
457
486
    else 
458
487
 
459
488
       ! Print-out whether transiesta is starting, or siesta is starting
460
 
       call ts_method_init( init_method == 2 )
 
489
       call ts_method_init( abs(init_method) == 3 )
461
490
 
462
491
    end if
463
492
 
559
588
    !                                 method used to read the DM/EDM
560
589
    !                                   0 == not read
561
590
    !                                   1 == .DM read
562
 
    !                                   2 == .TSDE read
 
591
    !                                   3 == .TSDE read
563
592
    ! *******************************************************************
564
593
 
565
594
    !> The spin-configuration that is used to determine the spin-order.
592
621
 
593
622
    ! The currently read stuff
594
623
    type(dSpData2D) :: DM_read
 
624
    type(Sparsity), pointer :: sp_read
595
625
    type(dSpData2D) :: EDM_read
596
626
 
597
627
    ! Signal the file has not been found.
614
644
 
615
645
       if ( TSDE_found ) then
616
646
          ! Signal we have read TSDE
617
 
          init_method = 2
 
647
          init_method = 3
618
648
 
619
649
          DM_found = .true.
620
650
 
670
700
    ! Density matrix size checks
671
701
    if ( DM_found ) then
672
702
 
 
703
      ! Specify the *correct* init_method
 
704
      ! In cases where the sparsity pattern changes, we assume this has to do with
 
705
      ! MD simulations and thus a geometry change.
 
706
      ! By default we do not allow mixing the first step when re-using a DM from another
 
707
      ! geometry. This is because of idempotency.
 
708
      sp_read => spar(DM_read)
 
709
      if ( .not. equivalent(sp, sp_read) ) then
 
710
        init_method = -init_method
 
711
      end if
 
712
 
673
713
      corrected_nsc = .false.
674
714
      if ( nsc_read(1) /= 0 .and. any(nsc /= nsc_read) ) then
675
715
 
959
999
!$OMP parallel default(shared), private(i1,i2)
960
1000
 
961
1001
      ! Initialize to 0
962
 
!$OMP do collapse(2)
963
1002
      do i2 = 1 , spin%DM
 
1003
!$OMP do
964
1004
         do i1 = 1 , nnz
965
1005
            DM(i1,i2) = 0._dp
966
1006
         end do
 
1007
!$OMP end do nowait
967
1008
      end do
968
 
!$OMP end do
969
 
 
 
1009
       
 
1010
!$OMP barrier
970
1011
      
971
1012
      ! Automatic, for non magnetic (nspin=1) or for Ferro or Antiferro
972
1013
!$OMP single
1164
1205
!$OMP parallel default(shared), private(i,ind,is,jo,gio,ia,io)
1165
1206
 
1166
1207
       ! Initialize to 0
1167
 
!$OMP do collapse(2)
1168
1208
       do is = 1 , spin%DM
 
1209
!$OMP do
1169
1210
          do i = 1 , nnz
1170
1211
             DM(i,is) = 0._dp
1171
1212
          end do
 
1213
!$OMP end do nowait
1172
1214
       end do
1173
 
!$OMP end do
1174
1215
 
 
1216
!$OMP barrier
 
1217
       
1175
1218
       ! Initialize the paramagnetic case
1176
1219
!$OMP single
1177
1220
       do io = 1 , no_l