Changeset 702


Ignore:
Timestamp:
04/26/12 10:24:56 (14 months ago)
Author:
micael
Message:
  • Removed the PPScheme variable and changed the PPComponents block so that the scheme is specified for each component. This means it is now possible to use different schemes at the same time.
  • Valence states higher in energy than the ones specified in the PPComponents block are now included properly.
  • Now only semi-core component should be specified for the MRPP scheme, as the code will find the corresponding valence state automatically. Also, the same cut-off radius will be used for both, as I realized that in fact it made no sense for it to be different.
Location:
trunk/src
Files:
5 edited

Legend:

Unmodified
Added
Removed

  • trunk/src/atom_ps_inc.F90

    r683 r702  
    3131      real(R8) :: j 
    3232      real(R8) :: rc 
     33      integer  :: scheme 
    3334    end type input_data 
    3435 
    3536    type channel 
    36       integer :: scheme 
    3737      type(states_batch_t) :: fold 
    38       real(R8),   pointer  :: rc(:) 
    39       type(qn_t), pointer  :: qn(:) 
     38      integer    :: scheme 
     39      real(R8)   :: rc 
     40      type(qn_t) :: qn 
    4041    end type channel 
    4142 
    42     integer :: i, j, k, n, scheme, ierr, n_cols, n_lines, n_states, n_channels, fold_size, nlcc, file_format 
     43    integer :: i, j, k, n, ierr, n_cols, n_lines, n_channels, nlcc, file_format 
    4344    integer(POINTER_SIZE) :: blk 
    4445    real(R8) :: z_val, tol, min_ev, max_ev, ev, norm, slope, rmatch, nlcc_rc 
    4546    character(3) :: unbound_trick 
    4647    character(10) :: label 
    47     character(40) :: scheme_name 
    48     type(qn_t) :: fold_qn 
    4948    real(R8), allocatable :: rc(:), core_density(:,:), ps_density(:,:), ps_v(:,:), vh(:,:), vxc(:,:), vxctau(:,:) 
    5049    type(input_data), allocatable :: id(:) 
     50    type(qn_t), allocatable :: qn(:) 
     51    type(state_t), pointer :: ae_state, state 
    5152    type(channel), allocatable :: channels(:) 
    52     type(qn_t), allocatable :: qn(:) 
    53     type(state_t), pointer :: ae_state, state, state2 
    54     type(states_batch_t), allocatable :: folds(:) 
    5553 
    5654    call push_sub("atom_create_ps") 
    … …  
    8179    ps_atm%integrator_dp = ae_atm%integrator_dp 
    8280 
     81 
     82                    !---Read input options---! 
     83 
    8384    !Print some information 
    8485    message(1) = "" 
    … …  
    8788    call write_info(2,unit=info_unit("pp")) 
    8889 
    89     !What scheme are we using to generate the PP? 
    90     call oct_parse_int('PPScheme', HAM, scheme) 
    91     select case (scheme) 
    92     case (HAM,TM,MRPP,MTM) 
    93     case (RTM) 
    94       if (ps_atm%wave_eq /= DIRAC) then 
    95         message(1) = "All electron relativistic calculation is needed in order to use" 
    96         message(2) = "the relativistic extension to the Troullier-Martins scheme ." 
    97         call write_fatal(2) 
    98       end if 
    99     case (RMRPP) 
    100       if (ps_atm%wave_eq /= DIRAC) then 
    101         message(1) = "All electron relativistic calculation is needed in order." 
    102         message(2) = "to use the relativistic extension to the MRPP scheme" 
    103         call write_fatal(2) 
    104       end if 
    105     case default 
    106       message(1) = "Illegal pseusopotential generation scheme." 
    107       call write_fatal(1) 
    108     end select 
    109  
    11090    !Get PP generation tolerance 
    11191    call oct_parse_float('PPCalculationTolerance', 1e-6_r8, tol) 
    … …  
    128108 
    129109      !Check for errors in the format 
    130       if (n_cols < 3 .or. n_cols > 4 .or. & 
    131            (ps_atm%wave_eq /= DIRAC .and. n_cols == 4)) then 
     110      if (n_cols < 4 .or. n_cols > 5 .or. & 
     111           (ps_atm%wave_eq /= DIRAC .and. n_cols == 5)) then 
    132112        write(message(1),'("There is something wrong at line ",I2," of the PPComponents block.")') i 
    133113        call write_fatal(1) 
    … …  
    140120      call oct_parse_block_int(blk, i-1, 1, id(i)%l) 
    141121 
    142       !Get cutoff radius (it is always the last column) 
    143       call oct_parse_block_double(blk, i-1, n_cols-1, id(i)%rc) 
    144       id(i)%rc = id(i)%rc*units_in%length%factor 
    145        
    146       if (ps_atm%wave_eq == DIRAC .and. n_cols == 4) then 
     122      if (ps_atm%wave_eq == DIRAC .and. n_cols == 5) then 
    147123        !The third column is the total angular momentum quantum number 
    148124        call oct_parse_block_double(blk, i-1, 2, id(i)%j) 
    … …  
    150126        id(i)%j = M_ZERO 
    151127      end if 
     128 
     129      !Get cutoff radius (it is always the column before the last) 
     130      call oct_parse_block_double(blk, i-1, n_cols-2, id(i)%rc) 
     131      id(i)%rc = id(i)%rc*units_in%length%factor 
     132 
     133      !Get scheme (it is always the last column) 
     134      call oct_parse_block_int(blk, i-1, n_cols-1, id(i)%scheme) 
    152135 
    153136    end do 
    … …  
    163146      do j = i+1, n_lines 
    164147        if (id(i)%n == id(j)%n .and. id(i)%l == id(j)%l .and. id(i)%j == id(j)%j) then 
    165           write(message(1),'("There is something wrong at line ",I2," of the PPComponents block.")') j 
     148          write(message(1),'("There is something wrong at line ",I2," of the PPComponents block.")') i 
    166149          message(2) = "There is at least another line with the same set of quantum numbers." 
    167150          call write_fatal(2) 
    168151        end if 
    169152      end do 
    170     end do 
    171  
    172     !Number of states used to generate the pseudopotentials and respective quantum numbers and cutoff radii 
    173     n_states = n_lines 
    174     if (ps_atm%wave_eq == DIRAC) n_states = n_states + count(id%l /= 0 .and. id%j == M_ZERO) 
    175  
    176     allocate(qn(n_states), rc(n_states)) 
     153      select case (id(i)%scheme) 
     154      case (HAM,TM,MRPP,MTM) 
     155      case (RTM,RMRPP) 
     156        if (ps_atm%wave_eq /= DIRAC) then 
     157          write(message(1),'("There is something wrong at line ",I2," of the PPComponents block.")') i 
     158          message(2) = "All electron relativistic calculation is needed in order to use" 
     159          message(3) = "the relativistic extensions to the TM or MRPP schemes." 
     160          call write_fatal(2) 
     161        end if 
     162      case default 
     163        write(message(1),'("There is something wrong at line ",I2," of the PPComponents block.")') i 
     164        message(2) = "Illegal pseusopotential generation scheme." 
     165        call write_fatal(2) 
     166      end select 
     167 
     168    end do 
     169 
     170 
     171                    !---Initialize data structures---! 
     172 
     173    !Allocate memory and copy some information 
     174    n_channels = n_lines 
     175    if (ps_atm%wave_eq == DIRAC) n_channels = n_channels + count(id%l /= 0 .and. id%j == M_ZERO) 
     176    allocate(channels(n_channels)) 
     177    do i = 1, n_channels 
     178      call states_batch_null(channels(i)%fold) 
     179    end do 
    177180    k = 0 
    178181    do i = 1, n_lines 
    179182      k = k + 1 
    180183      if (ps_atm%wave_eq == DIRAC .and. id(i)%j == M_ZERO) then 
    181         qn(k) = qn_init(id(i)%n, id(i)%l, M_ZERO, j=id(i)%l+M_HALF) 
    182         rc(k) = id(i)%rc 
     184        !Split "j=l +/- 1/2" states that were implicit in the input file 
     185        channels(k)%qn = qn_init(id(i)%n, id(i)%l, M_ZERO, j=id(i)%l+M_HALF) 
     186        channels(k)%rc = id(i)%rc 
     187        channels(k)%scheme = id(i)%scheme 
    183188        if (id(i)%l /= M_ZERO) then 
    184189          k = k + 1 
    185           qn(k) = qn_init(id(i)%n, id(i)%l, M_ZERO, j=id(i)%l-M_HALF) 
    186           rc(k) = id(i)%rc 
     190          channels(k)%qn = qn_init(id(i)%n, id(i)%l, M_ZERO, j=id(i)%l-M_HALF) 
     191          channels(k)%rc = id(i)%rc 
     192          channels(k)%scheme = id(i)%scheme 
    187193        end if 
    188194      else 
    189         qn(k) = qn_init(id(i)%n, id(i)%l, M_ZERO, j=id(i)%j) 
    190         rc(k) = id(i)%rc 
     195        channels(k)%qn = qn_init(id(i)%n, id(i)%l, M_ZERO, j=id(i)%j) 
     196        channels(k)%rc = id(i)%rc 
     197        channels(k)%scheme = id(i)%scheme 
    191198      end if 
    192199    end do 
    193200    deallocate(id) 
    194201 
    195     !Create states batch 
     202    !Create states batch and add states to the channels folds 
    196203    call states_batch_null(ps_atm%states) 
    197     do i = 1, n_states 
     204    do i = 1, n_channels 
    198205      do n = 1, states_batch_size(ae_atm%states) 
    199206        ae_state => states_batch_get(ae_atm%states, n) 
    200         if (state_qn(ae_state) == qn(i)) then 
     207        if (state_qn(ae_state) == channels(i)%qn) then 
    201208          allocate(state) 
    202209          call state_null(state) 
    203210          state = ae_state 
    204           call states_batch_add(ps_atm%states, state)           
     211          call states_batch_add(ps_atm%states, state) 
     212          call states_batch_add(channels(i)%fold, state) 
    205213          nullify(state) 
    206214          exit 
    … …  
    208216        if (n == states_batch_size(ae_atm%states)) then 
    209217          write(message(1),'("Orbital ",A," not found in all-electron calculation.")') & 
    210                trim(qn_label(qn(i))) 
     218               trim(qn_label(channels(i)%qn)) 
    211219          call write_fatal(1) 
    212220        end if 
    … …  
    214222    end do 
    215223 
     224    !Sanity check 
     225    if (n_channels /= states_batch_number_of_folds(ps_atm%states)) then 
     226      write(message(1),'("There is something wrong at the PPComponents block.")') 
     227      message(2) = "Only one orbital per angular momentum channel should be specified." 
     228    end if 
     229 
    216230    !Check core radius 
    217     do i = 1, states_batch_size(ps_atm%states) 
    218       if (rc(i) == M_ZERO) then 
    219         rc(i) = state_default_rc(states_batch_get(ps_atm%states, i), scheme) 
    220       end if 
    221       if (rc(i) == M_ZERO) then 
     231    do i = 1, n_channels 
     232      if (channels(i)%rc == M_ZERO) then 
     233        channels(i)%rc = state_default_rc(states_batch_get(channels(i)%fold, 1), channels(i)%scheme) 
     234      end if 
     235      if (channels(i)%rc == M_ZERO) then 
    222236        do j = 1, states_batch_size(ps_atm%states) 
    223           rc(i) = max(rc(i), state_default_rc(states_batch_get(ps_atm%states, j), scheme)) 
     237          channels(i)%rc = max(channels(i)%rc, state_default_rc(states_batch_get(ps_atm%states, j), channels(i)%scheme)) 
    224238        end do 
    225239      end if 
    … …  
    238252    end if 
    239253 
    240     !Write information 
    241     select case (scheme) 
    242     case (HAM) 
    243       scheme_name = "Hamann" 
    244     case (TM) 
    245       scheme_name = "Troullier-Martins" 
    246     case (RTM) 
    247       scheme_name = "Troullier-Martins Relativistic extension" 
    248     case (MRPP) 
    249       scheme_name = "Multireference Pseudopotentials (MRPP)" 
    250     case (RMRPP) 
    251       scheme_name = "MRPP Relativistic extension" 
    252     case (MTM) 
    253       scheme_name = "Modified Troullier-Martins" 
    254     end select 
    255     write(message(1),'(2X,"PP generation scheme : ",A)') trim(scheme_name) 
    256     write(message(2),'(2X,"Wavefunction info:")') 
    257     write(message(3),'(4X,"State   Occupation    Node radius   Peak radius   Default core radius")') 
    258     n = 3 
    259     do i = 1, states_batch_size(ps_atm%states) 
    260       state => states_batch_get(ps_atm%states, i) 
    261       n = n + 1 
    262       label = state_label(state) 
    263       select case (len_trim(label)) 
    264       case (2) 
    265         write(message(n),'(5X,A,5X,F7.2,7X,F8.3,6X,F8.3,10X,F8.3)') & 
    266              trim(label), state_charge(state), & 
    267              state_outermost_node(state), & 
    268              state_outermost_peak(state), & 
    269              state_default_rc(state, scheme) 
    270       case (5) 
    271         write(message(n),'(4X,A,3X,F7.2,7XF8.3,6X,F8.3,10X,F8.3)') & 
    272              trim(label), state_charge(state), & 
    273              state_outermost_node(state), & 
    274              state_outermost_peak(state), & 
    275              state_default_rc(state, scheme) 
    276       end select 
    277     end do 
    278     message(n+1) = "" 
    279     message(n+2) = "Pseudopotential Generation:" 
    280     call write_info(n+2) 
    281     call write_info(n+2,unit=info_unit("pp")) 
    282  
    283  
    284     !Divide valence states into channels 
    285     n_channels = states_batch_number_of_folds(ps_atm%states) 
    286     allocate(folds(n_channels), channels(n_channels)) 
     254    !We will include in the valence space all the states that are higher in  
     255    !energy than the states used to generate the pseudopotentials 
    287256    do i = 1, n_channels 
    288       call states_batch_null(folds(i)) 
    289       call states_batch_null(channels(i)%fold) 
    290     end do 
    291     call states_batch_split_folds(ps_atm%states, folds) 
    292  
    293     do i = 1, n_channels 
    294       !Sanity check 
    295       fold_size = states_batch_size(folds(i)) 
    296       if (scheme /= MRPP .and. fold_size > 1) then 
    297         write(message(1),'("There is something wrong at the PPComponents block.")') i 
    298         message(2) = "When not using the MRPP scheme, only one orbital can be used per" 
    299         message(3) = "angular momentum channel to generate the pseudopotentials." 
    300         call write_fatal(3) 
    301       else if (scheme == MRPP .and. fold_size > 2) then 
    302         write(message(1),'("There is something wrong at the PPComponents block.")') i 
    303         message(2) = "When using the MRPP scheme, at most two orbital can be used per" 
    304         message(3) = "angular momentum channel to generate the pseudopotentials." 
    305         call write_fatal(3) 
    306       end if 
    307  
    308       !Actual scheme used to generate the pseudopotential 
    309       if (scheme == MRPP .and. fold_size == 1) then 
    310         channels(i)%scheme = TM 
    311       elseif (scheme == RMRPP .and. fold_size == 1) then 
    312         channels(i)%scheme = RTM 
    313       else 
    314         channels(i)%scheme = scheme 
    315       end if 
    316  
    317       channels(i)%fold = folds(i) 
    318       call states_batch_sort(channels(i)%fold, SORT_EV) 
    319       allocate(channels(i)%qn(fold_size)) 
    320       allocate(channels(i)%rc(fold_size)) 
    321  
    322       k = 0 
    323       fold_qn = state_qn(states_batch_get(channels(i)%fold, 1)) 
    324       do j = 1, n_states 
    325         if (qn_equal_fold(fold_qn, qn(j))) then 
    326           k = k + 1 
    327           channels(i)%qn(k) = qn(j) 
    328           channels(i)%rc(k) = rc(j) 
    329         end if 
    330       end do 
    331  
    332       !We will include in the valence space all the states that are higher in  
    333       !energy than the states used to generate the pseudopotentials 
    334257      min_ev = maxval(states_batch_eigenvalues(channels(i)%fold)) 
    335258      do j = 1, states_batch_size(ae_atm%states) 
    336259        ae_state => states_batch_get(ae_atm%states, j) 
    337         if (any(channels(i)%qn == state_qn(ae_state))) cycle 
    338         if (state_eigenvalue(ae_state) > min_ev .and. qn_equal_fold(state_qn(ae_state), fold_qn)) then 
     260        if (channels(i)%qn == state_qn(ae_state)) cycle 
     261        if (state_eigenvalue(ae_state) > min_ev .and. qn_equal_fold(state_qn(ae_state), channels(i)%qn)) then 
    339262          allocate(state) 
    340263          call state_null(state) 
    … …  
    344267          nullify(state) 
    345268          exit 
    346         end if 
    347        
    348       end do 
    349  
    350     end do 
    351     deallocate(qn) 
     269        end if       
     270      end do 
     271 
     272      !Sanity check 
     273      if ( (channels(i)%scheme == MRPP .or.  channels(i)%scheme == RMRPP) .and. & 
     274           states_batch_size(channels(i)%fold) == 1) then 
     275        write(message(1),'("There is something wrong at the PPComponents block.")') 
     276        message(2) = "When using the MRPP scheme, at least two orbitals from the" 
     277        message(3) = "same angular momentum channel are needed in the valence space," 
     278        message(4) = "but only one was found." 
     279        call write_fatal(4) 
     280      end if 
     281 
     282    end do 
     283 
     284    !Write information 
     285    write(message(1),'(2X,"Wavefunction info:")') 
     286    write(message(2),'(4X,"State   Occupation    Node radius   Peak radius   Default core radius")') 
     287    n = 2 
     288    do j = 1, n_channels 
     289      call states_batch_sort(channels(j)%fold, SORT_EV) 
     290      do i = 1, states_batch_size(channels(j)%fold) 
     291        state => states_batch_get(channels(j)%fold, i) 
     292        n = n + 1 
     293        label = state_label(state) 
     294        select case (len_trim(label)) 
     295        case (2) 
     296          write(message(n),'(5X,A,5X,F7.2,7X,F8.3,6X,F8.3,10X,F8.3)') & 
     297               trim(label), state_charge(state), & 
     298               state_outermost_node(state), & 
     299               state_outermost_peak(state), & 
     300               state_default_rc(state, channels(j)%scheme) 
     301        case (5) 
     302          write(message(n),'(4X,A,3X,F7.2,7XF8.3,6X,F8.3,10X,F8.3)') & 
     303               trim(label), state_charge(state), & 
     304               state_outermost_node(state), & 
     305               state_outermost_peak(state), & 
     306               state_default_rc(state, channels(j)%scheme) 
     307        end select 
     308      end do 
     309    end do 
     310    message(n+1) = "" 
     311    message(n+2) = "Pseudopotential Generation:" 
     312    call write_info(n+2) 
     313    call write_info(n+2,unit=info_unit("pp")) 
    352314 
    353315    !Get core density 
    … …  
    362324    ps_v = M_ZERO     
    363325    do i = 1, n_channels 
    364       state  => states_batch_get(channels(i)%fold, 1) 
    365       if (channels(i)%scheme == MRPP .or. channels(i)%scheme == RMRPP) then 
    366         state2 => states_batch_get(channels(i)%fold, 2) 
    367         call state_psp_generation(ps_atm%m, channels(i)%scheme, ps_atm%wave_eq, tol, & 
    368              ae_atm%potential, ps_atm%integrator_sp, ps_atm%integrator_dp, ps_v(:, i), & 
    369              state, channels(i)%rc(1), state2, channels(i)%rc(2)) 
    370       else 
    371         call state_psp_generation(ps_atm%m, channels(i)%scheme, ps_atm%wave_eq, tol, & 
    372              ae_atm%potential, ps_atm%integrator_sp, ps_atm%integrator_dp, ps_v(:, i), & 
    373              state, channels(i)%rc(1)) 
    374       end if 
     326      call states_batch_psp_generation(channels(i)%fold, ps_atm%m, channels(i)%scheme, & 
     327           ps_atm%wave_eq, tol, ae_atm%potential, ps_atm%integrator_sp, & 
     328           ps_atm%integrator_dp, eigensolver, channels(i)%rc, ps_v(:,i)) 
    375329    end do 
    376330 
    … …  
    395349      if (nlcc == CC_FHI) then 
    396350        call oct_parse_float('CoreCorrectionCutoff', ps_atm%m%r(i), nlcc_rc) 
    397       endif 
     351      end if 
    398352     
    399353      !Initialize the nlcc part of xc_model 
    … …  
    435389 
    436390        if (channels(i)%scheme == HAM) then 
    437           rmatch = hamann_match_radius(ps_atm%m, channels(i)%rc(j)) 
     391          rmatch = hamann_match_radius(ps_atm%m, channels(i)%rc) 
    438392        else 
    439           rmatch = channels(i)%rc(j) 
     393          rmatch = channels(i)%rc 
    440394        end if 
    441395 
    … …  
    464418    !Start writing things to the pseudopotential file 
    465419    call oct_parse_int('PPOutputFileFormat', PSPIO_UPF, file_format) 
    466     call ps_io_init(ps_atm%m, ps_atm%wave_eq, scheme, ae_atm%z, & 
     420    call ps_io_init(ps_atm%m, ps_atm%wave_eq, channels(1)%scheme, ae_atm%z, & 
    467421         ps_atm%z, z_val, ps_atm%symbol, file_format) 
    468422    call potential_ps_io_set(ps_atm%potential) 
     423    allocate(rc(states_batch_size(ps_atm%states))) 
     424    do i = 1, states_batch_size(ps_atm%states) 
     425      state => states_batch_get(ps_atm%states, i) 
     426      do j = 1, n_channels 
     427        if (qn_equal_fold(channels(j)%qn, state_qn(state))) then 
     428          rc(i) = channels(j)%rc 
     429          exit 
     430        end if 
     431      end do 
     432    end do 
    469433    call states_batch_ps_io_set(ps_atm%states, ps_atm%m, rc) 
     434    deallocate(rc) 
    470435    call xc_ps_io_set(ps_atm%xc_model) 
     436 
     437    !Free memory 
     438    do i = 1, n_channels 
     439      call states_batch_end(channels(i)%fold) 
     440    end do 
     441    deallocate(channels) 
    471442 
    472443    call pop_sub() 

  • trunk/src/mrpp.F90

    r653 r702  
    4444 
    4545  subroutine mrpp_gen(qn0, qn1, ev0, ev1, wave_eq, rel, tol, m, ae_potential, & 
    46                       rc0, rc1, integrator, ps_v, ps_wf0, ps_wf0p, ps_wf1, ps_wf1p) 
     46                      rc, integrator, ps_v, ps_wf0, ps_wf0p, ps_wf1, ps_wf1p) 
    4747    !-----------------------------------------------------------------------! 
    4848    ! Generates the pseudo wave-functions and the corresponding             ! 
    … …  
    5858    !  m            - mesh                                                  ! 
    5959    !  ae_potential - all-electron potential                                ! 
    60     !  rc0          - semi-core state cutoff radius                         ! 
    61     !  rc1          - valence state cutoff radius                           ! 
     60    !  rc           - cutoff radius                                         ! 
    6261    !  integrator   - integrator object                                     ! 
    6362    !  ps_v         - pseudo-potential on the mesh                          ! 
    … …  
    7877    type(mesh_t),       intent(in)    :: m 
    7978    type(potential_t),  intent(in)    :: ae_potential 
    80     real(R8),           intent(in)    :: rc0, rc1 
     79    real(R8),           intent(in)    :: rc 
    8180    type(integrator_t), intent(inout) :: integrator 
    8281    real(R8),           intent(out)   :: ps_v(m%np) 
    … …  
    8786 
    8887    integer :: wf_dim, i 
    89     real(R8) :: wf0_rc0, wf0p_rc0, wf1_rc1, norm0 
     88    real(R8) :: wf0_rc, wf0p_rc, wf1_rc, norm0 
    9089    real(R8) :: rhs(9), c(9) 
    9190    real(R8), allocatable :: ae_wf0(:,:), ae_wf0p(:,:) 
    9291    real(R8), allocatable :: ae_wf1(:,:), ae_wf1p(:,:) 
    93     type(mesh_t) :: m0, m1 
     92    type(mesh_t) :: mrc 
    9493 
    9594    call push_sub("mrpp_gen") 
    9695 
    9796    ! 
    98     call mesh_null(m0) 
    99     call mesh_null(m1) 
     97    call mesh_null(mrc) 
    10098 
    10199    !Write core radii 
    102     write(message(1),'(4x,"Core radius:",1x,f7.3,5x,f7.3)') rc0, rc1 
     100    write(message(1),'(4x,"Core radius:",1x,f7.3)') rc 
    103101    call write_info(1,20) 
    104102    call write_info(1,unit=info_unit("pp")) 
    … …  
    113111    ae_wf1p = ps_wf1p 
    114112 
    115     !Semi-core wavefunctions and norm at rc0 
    116     wf0_rc0  = rc0*mesh_extrapolate(m, ae_wf0(:,1), rc0) 
    117     wf0p_rc0 = rc0*mesh_extrapolate(m, ae_wf0p(:,1), rc0) + wf0_rc0/rc0 
    118     if (wf0_rc0 < M_ZERO) then 
     113    !Semi-core wavefunctions and norm at rc 
     114    wf0_rc  = rc*mesh_extrapolate(m, ae_wf0(:,1), rc) 
     115    wf0p_rc = rc*mesh_extrapolate(m, ae_wf0p(:,1), rc) + wf0_rc/rc 
     116    if (wf0_rc < M_ZERO) then 
    119117      ae_wf0  = -ae_wf0 
    120118      ae_wf0p = -ae_wf0p 
    121       wf0_rc0  = -wf0_rc0 
    122       wf0p_rc0 = -wf0p_rc0 
     119      wf0_rc  = -wf0_rc 
     120      wf0p_rc = -wf0p_rc 
    123121    end if 
    124122    if (rel) then 
    125       norm0 = log(mesh_integrate(m, sum(ae_wf0**2,dim=2), b=rc0)) 
     123      norm0 = log(mesh_integrate(m, sum(ae_wf0**2,dim=2), b=rc)) 
    126124    else 
    127125      select case (wave_eq) 
    128126      case (SCHRODINGER, SCALAR_REL) 
    129         norm0 = log(mesh_integrate(m, ae_wf0(:,1)**2, b=rc0)) 
     127        norm0 = log(mesh_integrate(m, ae_wf0(:,1)**2, b=rc)) 
    130128      case (DIRAC) 
    131         norm0 = log(mesh_integrate(m, ae_wf0(:,1)**2, b=rc0) + mesh_integrate(m, ae_wf0(:,2)**2)) 
     129        norm0 = log(mesh_integrate(m, ae_wf0(:,1)**2, b=rc) + mesh_integrate(m, ae_wf0(:,2)**2)) 
    132130      end select 
    133131    end if 
    134132 
    135     !Valence wavefunctions at rc1 
    136     wf1_rc1  = rc1*mesh_extrapolate(m, ae_wf1(:,1), rc1) 
     133    !Valence wavefunctions at rc 
     134    wf1_rc  = rc*mesh_extrapolate(m, ae_wf1(:,1), rc) 
    137135 
    138136    !Get right hand side of equations to solve 
    139     call tm_equations_rhs(m, rc0, rel, qn0, ev0, wf0_rc0, wf0p_rc0, norm0, ae_potential, rhs(1:7)) 
    140     rhs(8) = wf1_rc1 
     137    call tm_equations_rhs(m, rc, rel, qn0, ev0, wf0_rc, wf0p_rc, norm0, ae_potential, rhs(1:7)) 
     138    rhs(8) = wf1_rc 
    141139    rhs(9) = M_ZERO 
    142140 
    143     !New meshes 
    144     m0 = m 
    145     m1 = m 
    146     call mesh_truncate(m0, rc0) 
    147     call mesh_truncate(m1, rc1) 
     141    !New mesh 
     142    mrc = m 
     143    call mesh_truncate(mrc, rc) 
    148144 
    149145    !Set the quantum numbers to their correct value 
    … …  
    152148 
    153149    !Solve the system of equations 
    154     call mrpp_solve_system(m, m0, m1, .false., wave_eq, integrator, ae_potential,& 
    155          rc0, rc1, qn0, qn1, ev0, ev1, rhs, tol, c, ps_wf0(1:m0%np,:), & 
    156          ps_wf0p(1:m0%np,:), ps_wf1(1:m1%np,:), ps_wf1p(1:m1%np,:), ps_v(1:m0%np)) 
     150    call mrpp_solve_system(m, mrc, .false., wave_eq, integrator, ae_potential,& 
     151         rc, qn0, qn1, ev0, ev1, rhs, tol, c, ps_wf0(1:mrc%np,:), & 
     152         ps_wf0p(1:mrc%np,:), ps_wf1(1:mrc%np,:), ps_wf1p(1:mrc%np,:), ps_v(1:mrc%np)) 
    157153 
    158154    !Write info 
    … …  
    168164 
    169165    !Compute the pseudo-wavefunction and the pseudopotential 
    170     ps_wf0(m0%np:m%np,1:wf_dim) = ae_wf0(m0%np:m%np,1:wf_dim) 
    171     ps_wf0p(m0%np:m%np,1:wf_dim) = ae_wf0p(m0%np:m%np,1:wf_dim) 
    172     ps_wf1(m1%np:m%np,1:wf_dim) = ae_wf1(m1%np:m%np,1:wf_dim) 
    173     ps_wf1p(m1%np:m%np,1:wf_dim) = ae_wf1p(m1%np:m%np,1:wf_dim) 
     166    ps_wf0(mrc%np:m%np,1:wf_dim) = ae_wf0(mrc%np:m%np,1:wf_dim) 
     167    ps_wf0p(mrc%np:m%np,1:wf_dim) = ae_wf0p(mrc%np:m%np,1:wf_dim) 
     168    ps_wf1(mrc%np:m%np,1:wf_dim) = ae_wf1(mrc%np:m%np,1:wf_dim) 
     169    ps_wf1p(mrc%np:m%np,1:wf_dim) = ae_wf1p(mrc%np:m%np,1:wf_dim) 
    174170    if (wave_eq == DIRAC .and. .not. rel) then 
    175       ps_wf0(m0%np:m%np,2) = M_ZERO 
    176       ps_wf0p(m0%np:m%np,2) = M_ZERO 
    177       ps_wf1(m1%np:m%np,2) = M_ZERO 
    178       ps_wf1p(m1%np:m%np,2) = M_ZERO 
     171      ps_wf0(mrc%np:m%np,2) = M_ZERO 
     172      ps_wf0p(mrc%np:m%np,2) = M_ZERO 
     173      ps_wf1(mrc%np:m%np,2) = M_ZERO 
     174      ps_wf1p(mrc%np:m%np,2) = M_ZERO 
    179175    end if 
    180     do i = m0%np, m%np 
     176    do i = mrc%np, m%np 
    181177      ps_v(i) = v(ae_potential, m%r(i), qn0) 
    182178    end do 
    … …  
    184180    !Free memory 
    185181    deallocate(ae_wf0, ae_wf0p, ae_wf1, ae_wf1p) 
    186     call mesh_end(m0) 
    187     call mesh_end(m1) 
     182    call mesh_end(mrc) 
    188183 
    189184    call pop_sub() 
    190185  end subroutine mrpp_gen 
    191186 
    192   subroutine mrpp_solve_system(m, m0, m1, rel, wave_eq, integrator, ae_potential,& 
    193                                rc0, rc1, qn0, qn1, ev0, ev1, rhs, tol, c, & 
     187  subroutine mrpp_solve_system(m, mrc, rel, wave_eq, integrator, ae_potential,& 
     188                               rc, qn0, qn1, ev0, ev1, rhs, tol, c, & 
    194189                               ps_wf0, ps_wf0p, ps_wf1, ps_wf1p, ps_v) 
    195190    !-----------------------------------------------------------------------! 
    … …  
    197192    !                                                                       ! 
    198193    !  m            - mesh                                                  ! 
    199     !  m0           - truncated mesh for semi-core state                    ! 
    200     !  m1           - truncated mesh for valence state                      ! 
     194    !  mrc          - truncated mesh                                        ! 
    201195    !  rel          - if true, use the relativistic extension to the scheme ! 
    202196    !  wave_eq      - wave equation to use                                  ! 
    203197    !  integrator   - integrator object                                     ! 
    204198    !  ae_potential - all-electron potential                                ! 
    205     !  rc0          - semi-core state cutoff radius                         ! 
    206     !  rc1          - valence state cutoff radius                           ! 
     199    !  rc           - cutoff radius                                         ! 
    207200    !  qn0          - quantum numbers of semi-core state                    ! 
    208201    !  qn1          - quantum numbers of valence state                      ! 
    … …  
    218211    !  ps_v         - pseudo-potential                                      ! 
    219212    !-----------------------------------------------------------------------! 
    220     type(mesh_t),       intent(in)    :: m, m0, m1 
     213    type(mesh_t),       intent(in)    :: m, mrc 
    221214    logical,            intent(in)    :: rel 
    222215    integer,            intent(in)    :: wave_eq 
    … …  
    224217    type(potential_t),  intent(in)    :: ae_potential 
    225218    type(qn_t),         intent(in)    :: qn0, qn1 
    226     real(R8),           intent(in)    :: rc0, rc1, ev0, ev1, rhs(9), tol 
     219    real(R8),           intent(in)    :: rc, ev0, ev1, rhs(9), tol 
    227220    real(R8),           intent(out)   :: c(9) 
    228     real(R8),           intent(inout) :: ps_wf0(:,:)  ! ps_wf0(m0%np, qn_wf_dim(qn)) 
    229     real(R8),           intent(inout) :: ps_wf0p(:,:) ! ps_wf0p(m0%np, qn_wf_dim(qn)) 
    230     real(R8),           intent(inout) :: ps_wf1(:,:)  ! ps_wf1(m1%np, qn_wf_dim(qn)) 
    231     real(R8),           intent(inout) :: ps_wf1p(:,:) ! ps_wf1p(m1%np, qn_wf_dim(qn)) 
    232     real(R8),           intent(out)   :: ps_v(m0%np) 
     221    real(R8),           intent(inout) :: ps_wf0(:,:)  ! ps_wf0(mrc%np, qn_wf_dim(qn)) 
     222    real(R8),           intent(inout) :: ps_wf0p(:,:) ! ps_wf0p(mrc%np, qn_wf_dim(qn)) 
     223    real(R8),           intent(inout) :: ps_wf1(:,:)  ! ps_wf1(mrc%np, qn_wf_dim(qn)) 
     224    real(R8),           intent(inout) :: ps_wf1p(:,:) ! ps_wf1p(mrc%np, qn_wf_dim(qn)) 
     225    real(R8),           intent(out)   :: ps_v(mrc%np) 
    233226 
    234227    real(R8) :: ldd, cc(4) 
    … …  
    239232    !Solve TM system: we will use it as a starting point 
    240233    c = M_ZERO 
    241     call tm_solve_system(m0, rc0, rel, qn0, ev0, rhs(1:7), tol, c(1:7), ps_wf0, ps_wf0p, ps_v, .true.) 
     234    call tm_solve_system(mrc, rc, rel, qn0, ev0, rhs(1:7), tol, c(1:7), ps_wf0, ps_wf0p, ps_v, .true.) 
    242235 
    243236    !Allocate memory 
    244     call mrpp_equations_init(m, m0, m1, rel, wave_eq, integrator, ae_potential, rc0, rc1, qn0, qn1, ev0, ev1, rhs) 
     237    call mrpp_equations_init(m, mrc, rel, wave_eq, integrator, ae_potential, rc, qn0, qn1, ev0, ev1, rhs) 
    245238 
    246239    !Solve system 

  • trunk/src/mrpp_equations.F90

    r653 r702  
    3636  logical            :: rel 
    3737  integer            :: wave_eq, wf_dim 
    38   real(R8)           :: rc0, rc1, ev0, ev1, rhs(9) 
    39   type(mesh_t)       :: m, m0, m1 
     38  real(R8)           :: rc, ev0, ev1, rhs(9) 
     39  type(mesh_t)       :: m, mrc 
    4040  type(potential_t)  :: ae_potential 
    4141  type(integrator_t) :: integrator 
    … …  
    6464contains 
    6565 
    66   subroutine mrpp_equations_init(mp, m0p, m1p, relp, wave_eqp, integratorp, & 
    67        ae_potentialp, rc0p, rc1p, qn0p, qn1p, ev0p, ev1p, rhsp) 
     66  subroutine mrpp_equations_init(mp, mrcp, relp, wave_eqp, integratorp, & 
     67       ae_potentialp, rcp, qn0p, qn1p, ev0p, ev1p, rhsp) 
    6868    !-----------------------------------------------------------------------! 
    6969    ! Initializes global data needed to solve the MRPP set of non-linear    ! 
    … …  
    7171    !                                                                       ! 
    7272    !  mp            - mesh                                                 ! 
    73     !  m0p           - truncated mesh for semi-core state                   ! 
    74     !  m1p           - truncated mesh for valence state                     ! 
     73    !  mrcp          - truncated mesh                                       ! 
    7574    !  relp          - if true, use the relativistic extension to the scheme! 
    7675    !  wave_eqp      - wave equation to use                                 ! 
    7776    !  integratorp   - integrator object                                    ! 
    7877    !  ae_potentialp - all-electron potential                               ! 
    79     !  rc0p          - semi-core state cutoff radius                        ! 
    80     !  rc1p          - valence state cutoff radius                          ! 
     78    !  rcp           - cutoff radius                                        ! 
    8179    !  qn0p          - quantum numbers of semi-core state                   ! 
    8280    !  qn1p          - quantum numbers of valence state                     ! 
    … …  
    8583    !  rhsp          - right-hand side of equations                         ! 
    8684    !-----------------------------------------------------------------------! 
    87     type(mesh_t),       intent(in) :: mp, m0p, m1p 
     85    type(mesh_t),       intent(in) :: mp, mrcp 
    8886    logical,            intent(in) :: relp 
    8987    integer,            intent(in) :: wave_eqp 
    … …  
    9189    type(potential_t),  intent(in) :: ae_potentialp 
    9290    type(qn_t),         intent(in) :: qn0p, qn1p 
    93     real(R8),           intent(in) :: rc0p, rc1p, ev0p, ev1p, rhsp(9) 
     91    real(R8),           intent(in) :: rcp, ev0p, ev1p, rhsp(9) 
    9492 
    9593    integer  :: i 
    … …  
    9997    ! 
    10098    call mesh_null(m) 
    101     call mesh_null(m0) 
    102     call mesh_null(m1) 
     99    call mesh_null(mrc) 
    103100    call potential_null(ae_potential) 
    104101 
    105102    ! 
    106     m = mp; m0 = m0p; m1 = m1p 
     103    m = mp; mrc = mrcp 
    107104    rel = relp 
    108105    wave_eq = wave_eqp 
    109106    integrator = integratorp 
    110107    ae_potential = ae_potentialp 
    111     rc0 = rc0p; rc1 = rc1p 
     108    rc = rcp 
    112109    qn0 = qn0p; qn1 = qn1p 
    113110    ev0 = ev0p; ev1 = ev1p 
    … …  
    115112 
    116113    !Initialize stuff to compute the polynomial 
    117     allocate(p0_r(m0%np, 9), p1_r(m0%np, 9), p2_r(m0%np, 9)) 
    118  
    119     p0_r(:, 1) = m0%r**2;  p1_r(:, 1) = m0%r;     p2_r(:, 1) = M_ONE; 
    120     p0_r(:, 2) = m0%r**4;  p1_r(:, 2) = m0%r**3;  p2_r(:, 2) = m0%r**2; 
    121     p0_r(:, 3) = m0%r**6;  p1_r(:, 3) = m0%r**5;  p2_r(:, 3) = m0%r**4; 
    122     p0_r(:, 4) = m0%r**8;  p1_r(:, 4) = m0%r**7;  p2_r(:, 4) = m0%r**6; 
    123     p0_r(:, 5) = m0%r**10; p1_r(:, 5) = m0%r**9;  p2_r(:, 5) = m0%r**8; 
    124     p0_r(:, 6) = m0%r**12; p1_r(:, 6) = m0%r**11; p2_r(:, 6) = m0%r**10; 
    125     p0_r(:, 7) = m0%r**14; p1_r(:, 7) = m0%r**13; p2_r(:, 7) = m0%r**12 
    126     p0_r(:, 8) = m0%r**16; p1_r(:, 8) = m0%r**15; p2_r(:, 8) = m0%r**14 
     114    allocate(p0_r(mrc%np, 9), p1_r(mrc%np, 9), p2_r(mrc%np, 9)) 
     115 
     116    p0_r(:, 1) = mrc%r**2;  p1_r(:, 1) = mrc%r;     p2_r(:, 1) = M_ONE; 
     117    p0_r(:, 2) = mrc%r**4;  p1_r(:, 2) = mrc%r**3;  p2_r(:, 2) = mrc%r**2; 
     118    p0_r(:, 3) = mrc%r**6;  p1_r(:, 3) = mrc%r**5;  p2_r(:, 3) = mrc%r**4; 
     119    p0_r(:, 4) = mrc%r**8;  p1_r(:, 4) = mrc%r**7;  p2_r(:, 4) = mrc%r**6; 
     120    p0_r(:, 5) = mrc%r**10; p1_r(:, 5) = mrc%r**9;  p2_r(:, 5) = mrc%r**8; 
     121    p0_r(:, 6) = mrc%r**12; p1_r(:, 6) = mrc%r**11; p2_r(:, 6) = mrc%r**10; 
     122    p0_r(:, 7) = mrc%r**14; p1_r(:, 7) = mrc%r**13; p2_r(:, 7) = mrc%r**12 
     123    p0_r(:, 8) = mrc%r**16; p1_r(:, 8) = mrc%r**15; p2_r(:, 8) = mrc%r**14 
    127124 
    128125    do i = 1, 8 
    … …  
    131128    end do 
    132129 
    133     p0_rc(1) = rc0**2;  p1_rc(1) = rc0;     p2_rc(1) = M_ONE;   p3_rc(1) = M_ZERO;  p4_rc(1) = M_ZERO 
    134     p0_rc(2) = rc0**4;  p1_rc(2) = rc0**3;  p2_rc(2) = rc0**2;  p3_rc(2) = rc0;     p4_rc(2) = M_ONE 
    135     p0_rc(3) = rc0**6;  p1_rc(3) = rc0**5;  p2_rc(3) = rc0**4;  p3_rc(3) = rc0**3;  p4_rc(3) = rc0**2 
    136     p0_rc(4) = rc0**8;  p1_rc(4) = rc0**7;  p2_rc(4) = rc0**6;  p3_rc(4) = rc0**5;  p4_rc(4) = rc0**4 
    137     p0_rc(5) = rc0**10; p1_rc(5) = rc0**9;  p2_rc(5) = rc0**8;  p3_rc(5) = rc0**7;  p4_rc(5) = rc0**6 
    138     p0_rc(6) = rc0**12; p1_rc(6) = rc0**11; p2_rc(6) = rc0**10; p3_rc(6) = rc0**9;  p4_rc(6) = rc0**8 
    139     p0_rc(7) = rc0**14; p1_rc(7) = rc0**13; p2_rc(7) = rc0**12; p3_rc(7) = rc0**11; p4_rc(7) = rc0**10 
    140     p0_rc(8) = rc0**16; p1_rc(8) = rc0**15; p2_rc(8) = rc0**14; p3_rc(8) = rc0**13; p4_rc(8) = rc0**12 
     130    p0_rc(1) = rc**2;  p1_rc(1) = rc;     p2_rc(1) = M_ONE;   p3_rc(1) = M_ZERO;  p4_rc(1) = M_ZERO 
     131    p0_rc(2) = rc**4;  p1_rc(2) = rc**3;  p2_rc(2) = rc**2;  p3_rc(2) = rc;     p4_rc(2) = M_ONE 
     132    p0_rc(3) = rc**6;  p1_rc(3) = rc**5;  p2_rc(3) = rc**4;  p3_rc(3) = rc**3;  p4_rc(3) = rc**2 
     133    p0_rc(4) = rc**8;  p1_rc(4) = rc**7;  p2_rc(4) = rc**6;  p3_rc(4) = rc**5;  p4_rc(4) = rc**4 
     134    p0_rc(5) = rc**10; p1_rc(5) = rc**9;  p2_rc(5) = rc**8;  p3_rc(5) = rc**7;  p4_rc(5) = rc**6 
     135    p0_rc(6) = rc**12; p1_rc(6) = rc**11; p2_rc(6) = rc**10; p3_rc(6) = rc**9;  p4_rc(6) = rc**8 
     136    p0_rc(7) = rc**14; p1_rc(7) = rc**13; p2_rc(7) = rc**12; p3_rc(7) = rc**11; p4_rc(7) = rc**10 
     137    p0_rc(8) = rc**16; p1_rc(8) = rc**15; p2_rc(8) = rc**14; p3_rc(8) = rc**13; p4_rc(8) = rc**12 
    141138 
    142139    p1_rc = p1_rc*FACTORS 
    … …  
    147144    ! 
    148145    wf_dim = qn_wf_dim(qn0) 
    149     allocate(ps_wf0(m0%np, wf_dim), ps_wf0p(m0%np, wf_dim)) 
    150     allocate(ps_wf1(m1%np, wf_dim), ps_wf1p(m1%np, wf_dim)) 
    151     allocate(ps_v(m0%np)) 
     146    allocate(ps_wf0(mrc%np, wf_dim), ps_wf0p(mrc%np, wf_dim)) 
     147    allocate(ps_wf1(mrc%np, wf_dim), ps_wf1p(mrc%np, wf_dim)) 
     148    allocate(ps_v(mrc%np)) 
    152149 
    153150  end subroutine mrpp_equations_init 
    … …  
    177174    f_x(1) = (xx(2)**2 + (M_TWO*qn0%l + M_FIVE)*xx(3) - rhs(6))/M_TEN 
    178175    if (rel) then 
    179       f_x(2) = log(mesh_integrate(m0, sum(ps_wf0**2, dim=2), b=rc0)) - rhs(7) 
     176      f_x(2) = log(mesh_integrate(mrc, sum(ps_wf0**2, dim=2), b=rc)) - rhs(7) 
    180177    else 
    181       f_x(2) = log(mesh_integrate(m0, ps_wf0(:,1)**2, b=rc0)) - rhs(7) 
     178      f_x(2) = log(mesh_integrate(mrc, ps_wf0(:,1)**2, b=rc)) - rhs(7) 
    182179    end if 
    183     f_x(3) = ps_wf1(m1%np,1)*rc1 - rhs(8) 
     180    f_x(3) = ps_wf1(mrc%np,1)*rc - rhs(8) 
    184181    f_x(4) = ldd - rhs(9) 
    185182 
    … …  
    214211    !-----------------------------------------------------------------------! 
    215212    real(R8), intent(in)  :: c(9) 
    216     real(R8), intent(out) :: wf0(m0%np, wf_dim), wf0p(m0%np, wf_dim) 
    217     real(R8), intent(out) :: pot(m0%np) 
    218     real(R8), intent(out) :: wf1(m1%np, wf_dim), wf1p(m1%np, wf_dim) 
     213    real(R8), intent(out) :: wf0(mrc%np, wf_dim), wf0p(mrc%np, wf_dim) 
     214    real(R8), intent(out) :: pot(mrc%np) 
     215    real(R8), intent(out) :: wf1(mrc%np, wf_dim), wf1p(mrc%np, wf_dim) 
    219216    real(R8), intent(out) :: ldd 
    220217 
    … …  
    226223 
    227224    !Compute the polynomial and its derivatives 
    228     allocate(p(m0%np), pp(m0%np), ppp(m0%np)) 
     225    allocate(p(mrc%np), pp(mrc%np), ppp(mrc%np)) 
    229226    p   = c(1) 
    230227    pp  = M_ZERO 
    … …  
    236233    end do 
    237234 
    238     pot = ev0 + (real(qn0%l,r8) + M_ONE)/m0%r*pp + (ppp + pp**2)/M_TWO 
    239     wf0(:,1) = m0%r**qn0%l*exp(p) 
     235    pot = ev0 + (real(qn0%l,r8) + M_ONE)/mrc%r*pp + (ppp + pp**2)/M_TWO 
     236    wf0(:,1) = mrc%r**qn0%l*exp(p) 
    240237    if (qn0%l == 0) then 
    241238      wf0p(:,1) = pp*exp(p) 
    242239    else 
    243       wf0p(:,1) = real(qn0%l,R8)*m0%r**(qn0%l-1)*exp(p) + m0%r**qn0%l*pp*exp(p) 
     240      wf0p(:,1) = real(qn0%l,R8)*mrc%r**(qn0%l-1)*exp(p) + mrc%r**qn0%l*pp*exp(p) 
    244241    end if 
    245242 
    … …  
    247244    if (rel) then 
    248245      !Relativistic correction to the pseudopotential 
    249       call potential_init(ps_potential, m0, pot) 
    250       do i = 1, m0%np 
    251         pot(i) = pot(i) + (pot(i) - ev0)**2/M_TWO/M_C2 + dvdr(ps_potential, m0%r(i), qn0)/M_FOUR/M_C2*(wf0p(i,1)/wf0(i,1)& 
    252                  + (qn0%k + M_ONE)/m0%r(i)) 
     246      call potential_init(ps_potential, mrc, pot) 
     247      do i = 1, mrc%np 
     248        pot(i) = pot(i) + (pot(i) - ev0)**2/M_TWO/M_C2 + dvdr(ps_potential, mrc%r(i), qn0)/M_FOUR/M_C2*(wf0p(i,1)/wf0(i,1)& 
     249                 + (qn0%k + M_ONE)/mrc%r(i)) 
    253250      end do 
    254251      call potential_end(ps_potential) 
    255252 
    256253      !Minor component 
    257       wf0(:,2) = M_C*((real(qn0%l,R8) + M_ONE + qn0%k)/m0%r + pp)*wf0(:,1)/(M_TWO*M_C2 - pot + ev0) 
     254      wf0(:,2) = M_C*((real(qn0%l,R8) + M_ONE + qn0%k)/mrc%r + pp)*wf0(:,1)/(M_TWO*M_C2 - pot + ev0) 
    258255      wf0p(:,2) = mesh_derivative(m, wf0(:,2)) 
    259256    else 
    … …  
    265262    ! 
    266263    allocate(vpp(m%np), wf(m%np, wf_dim), wfp(m%np, wf_dim)) 
    267     vpp(1:m0%np) = pot 
    268     do i = m0%np+1, m%np 
     264    vpp(1:mrc%np) = pot 
     265    do i = mrc%np+1, m%np 
    269266      vpp(i) = v(ae_potential, m%r(i), qn1) 
    270267    end do 
    … …  
    280277    end if 
    281278    do i = 1, wf_dim 
    282       call mesh_transfer(m, wf(:,i), m1, wf1(:,i)) 
    283       call mesh_transfer(m, wfp(:,i), m1, wf1p(:,i)) 
     279      call mesh_transfer(m, wf(:,i), mrc, wf1(:,i)) 
     280      call mesh_transfer(m, wfp(:,i), mrc, wf1p(:,i)) 
    284281    end do 
    285282    ldd = ldd*10.0 
    … …  
    315312    write(message(9),'(8X,"c2**2 + c4(2l+5)                  =",1x,ES16.9e2)') f(1) 
    316313    write(message(10),'(8X,"ln(norm_sc)[AE] - ln(norm_sc)[PP] =",1x,ES16.9e2)') f(2) 
    317     write(message(11),'(8X,"R(rc1)[AE] - R(rc1)[PP]           =",1x,ES16.9e2)') f(3) 
    318     write(message(12),'(8X,"R''/R(rc1)[AE] - R''/R(rc1)[PP]     =",1x,ES16.9e2)') f(4) 
     314    write(message(11),'(8X,"R(rc)[AE] - R(rc)[PP]             =",1x,ES16.9e2)') f(3) 
     315    write(message(12),'(8X,"R''/R(rc)[AE] - R''/R(rc)[PP]       =",1x,ES16.9e2)') f(4) 
    319316    message(13) = "" 
    320317    call write_info(13, 30) 
    … …  
    331328    deallocate(ps_wf0, ps_wf0p, ps_wf1, ps_wf1p, ps_v) 
    332329    call mesh_end(m) 
    333     call mesh_end(m0) 
    334     call mesh_end(m1) 
     330    call mesh_end(mrc) 
    335331    call potential_end(ae_potential) 
    336332 

  • trunk/src/states.F90

    r700 r702  
    962962  subroutine state_psp_generation(m, scheme, wave_eq, tol, ae_potential, & 
    963963                                  integrator_sp, integrator_dp, ps_v, state1, & 
    964                                   rc1, state2, rc2) 
     964                                  rc, state2) 
    965965    !-----------------------------------------------------------------------! 
    966966    ! Generate the pseudo wave-functions and the pseudo-potential using     ! 
    … …  
    976976    !  ps_v          - pseudo-potential on the mesh                         ! 
    977977    !  state1        - pseudo-state 1 (either valence or semi-core)         ! 
    978     !  rc1           - core radius of state 1                               ! 
     978    !  rc            - core radius                                          ! 
    979979    !  state2        - valence state when state 1 is a semi-core state      ! 
    980     !  rc2           - core radius of state 2                               ! 
    981980    !-----------------------------------------------------------------------! 
    982981    type(mesh_t),       intent(in)    :: m 
    … …  
    987986    real(R8),           intent(out)   :: ps_v(m%np) 
    988987    type(state_t),      intent(inout) :: state1 
    989     real(R8),           intent(in)    :: rc1 
     988    real(R8),           intent(in)    :: rc 
    990989    type(state_t),      intent(inout), optional :: state2 
    991     real(R8),           intent(in),    optional :: rc2 
    992  
    993     integer :: scheme_ 
     990 
    994991    character(len=10) :: label 
    995992 
    996993    call push_sub("state_psp_generation") 
    997994 
    998     scheme_ = scheme 
    999995    if (present(state2)) then 
    1000       ASSERT(present(rc2)) 
    1001     else 
    1002       select case (scheme_) 
    1003       case (MRPP) 
    1004         scheme_ = TM 
    1005       case (RMRPP) 
    1006         scheme_ = RTM 
    1007       end select 
     996      ASSERT(scheme == MRPP .or. scheme == RMRPP) 
    1008997    end if 
    1009998 
    … …  
    10181007    call write_info(1,unit=info_unit("pp")) 
    10191008 
    1020     select case (scheme_) 
     1009    select case (scheme) 
    10211010    case (HAM) !Hamann scheme 
    10221011      call hamann_gen(state1%qn, state1%ev, wave_eq, tol, m, ae_potential, & 
    1023                       rc1, integrator_sp, integrator_dp, ps_v, & 
     1012                      rc, integrator_sp, integrator_dp, ps_v, & 
    10241013                      state1%wf, state1%wfp) 
    10251014 
    10261015    case (TM) !Troullier-Martins scheme 
    1027       if (present(state2)) then 
    1028         !Get valence wave-function first 
    1029         call tm_gen(state2%qn, state2%ev, wave_eq, .false., tol, m, & 
    1030                     ae_potential, rc2, ps_v, state2%wf, state2%wfp) 
    1031         state2%qn%n = state2%qn%n + 1 
    1032       end if 
    1033  
    10341016      call tm_gen(state1%qn, state1%ev, wave_eq, .false., tol, m, & 
    1035                   ae_potential, rc1, ps_v, state1%wf, state1%wfp) 
     1017                  ae_potential, rc, ps_v, state1%wf, state1%wfp) 
    10361018       
    10371019    case (RTM) !Relativistic extension of the Troullier-Martins scheme 
    1038       if (present(state2)) then 
    1039         !Get valence wave-function first 
    1040         call tm_gen(state2%qn, state2%ev, wave_eq, .true., tol, m, & 
    1041                     ae_potential, rc2, ps_v, state2%wf, state2%wfp) 
    1042         state2%qn%n = state2%qn%n + 1 
    1043       end if 
    1044  
    10451020      call tm_gen(state1%qn, state1%ev, wave_eq, .true., tol, m, & 
    1046                   ae_potential, rc1, ps_v, state1%wf, state1%wfp) 
     1021                  ae_potential, rc, ps_v, state1%wf, state1%wfp) 
    10471022 
    10481023    case (MRPP) !Multireference Pseudopotentials 
    10491024      call mrpp_gen(state1%qn, state2%qn, state1%ev, state2%ev, & 
    1050                     wave_eq, .false., tol, m, ae_potential, rc1, & 
    1051                     rc2, integrator_dp, ps_v, state1%wf, state1%wfp, & 
     1025                    wave_eq, .false., tol, m, ae_potential, rc, & 
     1026                    integrator_dp, ps_v, state1%wf, state1%wfp, & 
    10521027                    state2%wf, state2%wfp) 
    10531028 
    10541029    case (RMRPP) !Multireference Pseudopotentials 
    10551030      call mrpp_gen(state1%qn, state2%qn, state1%ev, state2%ev, & 
    1056                     wave_eq, .true., tol, m, ae_potential, rc1, & 
    1057                     rc2, integrator_dp, ps_v, state1%wf, state1%wfp, & 
     1031                    wave_eq, .true., tol, m, ae_potential, rc, & 
     1032                    integrator_dp, ps_v, state1%wf, state1%wfp, & 
    10581033                    state2%wf, state2%wfp) 
    10591034 

  • trunk/src/states_batch.F90

    r699 r702  
    9797            states_batch_sort, & 
    9898            states_batch_smearing, & 
     99            states_batch_psp_generation, & 
    99100            states_batch_ld_test, & 
    100101            states_batch_output_configuration, & 
    … …  
    970971    call pop_sub() 
    971972  end subroutine states_batch_smearing 
     973 
     974  subroutine states_batch_psp_generation(batch, m, scheme, wave_eq, tol, ae_potential, & 
     975       integrator_sp, integrator_dp, eigensolver, rc, ps_v) 
     976    !-----------------------------------------------------------------------! 
     977    !-----------------------------------------------------------------------! 
     978    type(states_batch_t), intent(inout) :: batch 
     979    type(mesh_t),         intent(in)    :: m 
     980    integer,              intent(in)    :: scheme, wave_eq 
     981    real(R8),             intent(in)    :: tol 
     982    type(potential_t),    intent(in)    :: ae_potential 
     983    type(integrator_t),   intent(inout) :: integrator_sp, integrator_dp 
     984    type(eigensolver_t),  intent(in)    :: eigensolver 
     985    real(R8),             intent(in)    :: rc 
     986    real(R8),             intent(out)   :: ps_v(m%np) 
     987 
     988    integer :: n, i 
     989    type(state_t), pointer :: state, state2 
     990    type(potential_t) :: ps_potential 
     991    type(states_batch_t) :: other_states 
     992 
     993    call push_sub("states_batch_psp_generation") 
     994 
     995    state => states_batch_get(batch, 1) 
     996    if (scheme == MRPP .or. scheme == RMRPP) then 
     997      n = 2 
     998      state2 => states_batch_get(batch, 2) 
     999      call state_psp_generation(m, scheme, wave_eq, tol, ae_potential, & 
     1000           integrator_sp, integrator_dp, ps_v, state, rc, state2) 
     1001    else 
     1002      n = 1 
     1003      call state_psp_generation(m, scheme, wave_eq, tol, ae_potential, & 
     1004           integrator_sp, integrator_dp, ps_v, state, rc) 
     1005    end if 
     1006     
     1007 
     1008    !We need to solve the wave-equation with the pseudopotential for the other states in the batch 
     1009    if (n < states_batch_size(batch)) then 
     1010      call states_batch_null(other_states) 
     1011      call potential_null(ps_potential) 
     1012      call potential_init(ps_potential, m, ps_v) 
     1013 
     1014      do i = n+1, states_batch_size(batch) 
     1015        call states_batch_add(other_states, states_batch_get(batch, i)) 
     1016      end do 
     1017 
     1018      call states_batch_eigensolve(other_states, m, wave_eq, ps_potential, integrator_dp, integrator_sp, eigensolver) 
     1019 
     1020      call potential_end(ps_potential) 
     1021      call states_batch_end(other_states) 
     1022    end if 
     1023 
     1024    call pop_sub() 
     1025  end subroutine states_batch_psp_generation 
    9721026 
    9731027  subroutine states_batch_ld_test(batch, m, ae_potential, ps_potential, integrator, r, de, emin, emax) 
Note: See TracChangeset for help on using the changeset viewer.