← Back to branch summary

~jose-soler/siesta/newdynamics

« back to all changes in this revision

Viewing changes to Src/m_ts_electype.F90

  • Committer: Nick Papior
  • Date: 2018-10-18 20:36:17 UTC
  • mfrom: (560.1.457 4.1)
  • Revision ID: nickpapior@gmail.com-20181018203617-e2gk7g5vko9dmu55
Merged r1017, RE-SE fix for user-handling Hartree fix

Show diffs side-by-side

added added

removed removed

Lines of Context:
Src/m_ts_electype.F90
91
91
     ! chemical potential of the electrode
92
92
     type(ts_mu), pointer :: mu => null()
93
93
     ! infinity direction
94
 
     integer :: inf_dir = INF_NEGATIVE
 
94
     integer :: inf_dir = INF_POSITIVE
95
95
     ! transport direction (determines H01)
96
96
     ! And is considered with respect to the electrode direction...
97
97
     !  t_dir is with respect to the electrode unit-cell
399
399
           if ( len_trim(ln) /= 1 ) then
400
400
             call die(trim(tmp))
401
401
           end if
402
 
 
403
402
           ln = trim(ln) // fdf_bnames(pline,3)
404
 
 
405
403
         end if
406
404
 
407
 
         ! preset
408
 
         this%t_dir = 0           
 
405
         ! preset for checks of arguments
 
406
         this%t_dir = 0
 
407
         i = 0
409
408
 
410
409
         ! now for testing
411
410
         if ( ln(1:1) == '+' ) then
 
411
           i = 1
412
412
           this%inf_dir = INF_POSITIVE
 
413
           ln = ln(2:)
413
414
         else if ( ln(1:1) == '-' ) then
 
415
           i = 1
414
416
           this%inf_dir = INF_NEGATIVE
415
 
           ! The 3 below cases are *special* in the sense that they require the user
416
 
           ! to supply the GF files (TranSiesta/TBtrans cannot calculate the self-energies of
417
 
           ! real-space Green functions).
418
 
         else if ( leqi(ln,'ab') .or. leqi(ln, 'ba') ) then
 
417
           ln = ln(2:)
 
418
         end if
 
419
           
 
420
         ! The 3 below cases are *special* in the sense that they require the user
 
421
         ! to supply the GF files (TranSiesta/TBtrans cannot calculate the self-energies of
 
422
         ! real-space Green functions).
 
423
         if ( leqi(ln,'ab') .or. leqi(ln, 'ba') ) then
419
424
           this%t_dir = 6 ! Voigt notation
420
425
         else if ( leqi(ln,'ac') .or. leqi(ln, 'ca') ) then
421
426
           this%t_dir = 5
422
427
         else if ( leqi(ln,'bc') .or. leqi(ln, 'cb') ) then
423
428
           this%t_dir = 4
 
429
         else if ( leqi(ln,'c') .or. leqi(ln,'a3') ) then
 
430
           this%t_dir = 3
 
431
         else if ( leqi(ln,'b') .or. leqi(ln,'a2') ) then
 
432
           this%t_dir = 2
 
433
         else if ( leqi(ln,'a') .or. leqi(ln,'a1') ) then
 
434
           this%t_dir = 1
424
435
         else
 
436
           ! Simply a wrong argument (lattice vector)
425
437
           call die(trim(tmp))
426
438
         end if
427
439
 
428
 
         if ( this%t_dir > 3 ) then
429
 
           ! This flag has no influence
430
 
           this%inf_dir = INF_NEGATIVE
431
 
         else
432
 
           ! figure out direction...
433
 
           ln = ln(2:)
434
 
           if ( leqi(ln,'a') .or. leqi(ln,'a1') ) then
435
 
             this%t_dir = 1
436
 
           else if ( leqi(ln,'b') .or. leqi(ln,'a2') ) then
437
 
             this%t_dir = 2
438
 
           else if ( leqi(ln,'c') .or. leqi(ln,'a3') ) then
439
 
             this%t_dir = 3
440
 
           else
441
 
             call die(trim(tmp))
442
 
           end if
 
440
         ! In case only a single lattice vector has been specified, in that
 
441
         ! case the sign of the semi-infinite direction *must* be used!
 
442
         if ( this%t_dir <= 3 .and. i == 0 ) then
 
443
           call die(trim(tmp))
443
444
         end if
444
445
         info(2) = .true.
445
446
          
879
880
    ! Calculate planes of the electrodes
880
881
    select case ( this%t_dir )
881
882
    case ( 4 ) ! B-C
882
 
      p = this%cell(:,2)
883
 
      p = this%cell(:,3)
 
883
      call cross(this%cell(:,1), this%cell(:,2), p)
 
884
      contrib = VNORM(p)
 
885
      call cross(this%cell(:,1), this%cell(:,3), p)
 
886
      if ( contrib > VNORM(p) ) then
 
887
        ! the area on A-B is biggest, hence the normal plane is along C
 
888
        p = this%cell(:,3)
 
889
      else
 
890
        p = this%cell(:,2)
 
891
      end if
884
892
    case ( 5 ) ! A-C
885
 
      p = this%cell(:,1)
886
 
      p = this%cell(:,3)
 
893
      call cross(this%cell(:,2), this%cell(:,1), p)
 
894
      contrib = VNORM(p)
 
895
      call cross(this%cell(:,2), this%cell(:,3), p)
 
896
      if ( contrib > VNORM(p) ) then
 
897
        p = this%cell(:,3)
 
898
      else
 
899
        p = this%cell(:,1)
 
900
      end if
887
901
    case ( 6 ) ! A-B
888
 
      p = this%cell(:,1)
889
 
      p = this%cell(:,2)
 
902
      call cross(this%cell(:,3), this%cell(:,1), p)
 
903
      contrib = VNORM(p)
 
904
      call cross(this%cell(:,3), this%cell(:,2), p)
 
905
      if ( contrib > VNORM(p) ) then
 
906
        p = this%cell(:,2)
 
907
      else
 
908
        p = this%cell(:,1)
 
909
      end if
890
910
    case default
891
911
      p = this%cell(:,this%t_dir)
892
912
    end select
1552
1572
      ! Calculate pqosition in the grid
1553
1573
      call dgesv(3,1,LHS,3,idx,RHS,3,i)
1554
1574
      if ( i /= 0 ) then
1555
 
         call die('ts_voltage: Could not solve linear system.')
 
1575
         call die('ts_electype: Could not solve linear system.')
1556
1576
      end if
1557
1577
      
1558
1578
      ! Convert to integer position