- it is possible to specify j=l-1/2 and j=l+1/2 states for relativistic

calculations (rel=2) using "config" to pass the electronic configuratio
- more checks on what is read in input
- some cleanup


git-svn-id: http://qeforge.qe-forge.org/svn/q-e/trunk/espresso@2584 c92efa57-630b-4861-b058-cf58834340f0

atomic/el_config.f90

@@ -146,8 +146,10 @@ do n=nwfc+1,nwf
       if (oc(n) > 2*(2*ll(n)+1)) &
    &     call errore('el_config','wrong occupancy:'//occup,n)
 enddo
+!
 ! pseudopotentials: N quantum number .ne. wavefunction label
 ! we find the lowest N for each L and reassign N as it
+!
 if (.not.all_elec) then
    do l=0,3
       n0=1000

atomic/ld1_readin.f90

@@ -169,11 +169,11 @@ subroutine ld1_readin
      nspin = 1
   else if(lsd == 1) then
      nspin = 2
+     if (rel == 2) call errore('ld1_readin', &
+       &    'local spin density and spin-orbit not allowed',1)
   else
      call errore('ld1_readin','lsd not correct',1)
   endif
-  if (rel == 2 .and. lsd == 1) call errore('ld1_readin', &
-       &    'local spin density and spin-orbit not allowed',1)
 
   if (config == ' ') then
      call read_config (rel, lsd, nwf, el, nn, ll, oc, isw, jj)
@@ -182,12 +182,11 @@ subroutine ld1_readin
      !
      ! check same labels corresponding to different spin or j value
      !
-     do n=1,nwf  
+     jj(1:nwf)=0.d0
+     do n=1,nwf 
         do i=n+1,nwf  
            if (el(i) == el(n)) then
-              if ( lsd==0 ) then
-                 call errore('ld1_readin',el(i)//' appears twice',i)
-              else if (rel == 2) then
+              if (rel == 2) then
                  if (ll(n) > 0) then
                     jj(n) = ll(n) + (isw(n)-1.5)
                     jj(i) = ll(i) + (isw(i)-1.5)
@@ -198,12 +197,23 @@ subroutine ld1_readin
                  else
                     call errore('ld1_readin',el(i)//' appears twice',i)
                  end if
+              else if ( lsd==0 ) then
+                 call errore('ld1_readin',el(i)//' appears twice',i)
               end if
            end if
         end do
      end do
+     if (rel == 2) isw(1:nwf)=1 
   end if
   !
+  !  In the spin polarized or relativistic case adjust the occupations
+  !
+  if (lsd == 1) then
+     call occ_spin(nwf,nwfx,el,nn,ll,oc,isw)
+  else if (rel == 2) then
+     call occ_spinorb(nwf,nwfx,el,nn,ll,jj,oc,isw)
+  endif
+  !
   ! generate the radial grid - note that if iswitch = 2 the radial grid
   ! is not generated but read from the pseudopotential file
   !
@@ -212,19 +222,12 @@ subroutine ld1_readin
      rhoc=0.0_dp
   endif
   !
-  !  In the spin polarized case adjust the occupations
-  !
-  if (lsd == 1)  call occ_spin(nwf,nwfx,el,nn,ll,oc,isw)
-  !
-  !  If vdW calculation, read the input parameters
-  !
   if (iswitch == 1) then
      !
      !    no more data needed for AE calculations
      !
      return
-     !
-     
+     !     
   else if (iswitch == 3) then
      !
      !    reading input for PP generation
@@ -243,13 +246,16 @@ subroutine ld1_readin
 500  call errore('ld1_readin','reading inputp',abs(ios))
 
      if (lloc < 0 .and. rcloc <=0.0_dp) &
-          call errore('ld1_readin','rcloc is negative',1)
+          call errore('ld1_readin','rcloc must be positive',1)
      if (pseudotype < 1.or.pseudotype > 3) &
           call errore('ld1_readin','specify correct pseudotype',1)
      !
      call read_psconfig (rel, lsd, nwfs, els, nns, lls, ocs, &
           isws, jjs, enls, rcut, rcutus )
      !
+     if (rel==2) call occ_spinorbps &
+          (nwfs,nwfsx,els,nns,lls,jjs,ocs,rcut,rcutus,enls,isws)
+     !
      lmax = maxval(lls(1:nwfs))
      !
      zdum = zed
@@ -303,10 +309,6 @@ subroutine ld1_readin
   jjts=0.0_dp
   jjtsc=0.0_dp
   
-  do n=1,nwf
-     oc_old(n)=oc(n)
-  enddo
-
   nconf=1
   configts=' '
 
@@ -369,6 +371,11 @@ subroutine ld1_readin
      if (lsd == 1) then
         call occ_spin(nwftsc(nc),nwfsx,eltsc(1,nc),nntsc(1,nc),lltsc(1,nc),&
              octsc(1,nc), iswtsc(1,nc)) 
+     else if (rel == 2) then
+        call occ_spinorb(nwftsc(nc),nwfsx,eltsc(1,nc), &
+             &  nntsc(1,nc),lltsc(1,nc),jjtsc(1,nc),octsc(1,nc),iswtsc(1,nc))
+     else
+        jjtsc=0.0_dp
      endif
   end do
   !
@@ -542,6 +549,7 @@ subroutine read_config(rel, lsd, nwf, el, nn, ll, oc, isw, jj)
   !
   do n=1,nwf  
      if (rel < 2) then
+        jj(n) = 0.0_dp
         if (lsd == 0) then
            read(5,*,err=20,end=20,iostat=ios) &
                 el(n), nn(n), ll(n), oc(n)
@@ -559,12 +567,25 @@ subroutine read_config(rel, lsd, nwf, el, nn, ll, oc, isw, jj)
              el(n), nn(n), ll(n), oc(n), jj(n)
         isw(n)=1
         if ((abs(ll(n)+0.5_dp-jj(n)) > 1.e-3_dp) .and. &
-            (abs(ll(n)-0.5_dp-jj(n)) > 1.e-3_dp) .and. abs(jj(n)) > 1.e-3_dp)  &
+            (abs(ll(n)-0.5_dp-jj(n)) > 1.e-3_dp) .and. abs(jj(n)) > 1.e-3_dp) &
             call errore('read_config','jj wrong',n)
         if (oc(n) > (2.0_dp*jj(n)+1.0_dp) .and. abs(jj(n)) > 1e-3_dp) &
              call errore('read_config','occupations wrong',n)
 22      call errore('read_config','reading orbital (rel)',abs(ios))
      endif
+     !
+     ! Check: no two same wavefunctions
+     !
+     do ncheck=1,n-1
+        if ( el(ncheck) == el(n) .and. isw(ncheck) == isw(n) .and. &
+             jj(ncheck) == jj(n) ) then
+           call errore('read_config', &
+                'same wavefunction '//el(n)//' appears twice',n)
+        endif
+     enddo
+     !
+     ! More sanity checks
+     !
      write(label,'(a2)') el(n)
      read (label,'(i1)') ncheck
      if (ncheck /= nn(n)  .or. &

atomic/ld1_setup.f90

@@ -18,7 +18,7 @@ subroutine ld1_setup
   use funct, only : get_iexch, dft_is_meta !, set_dft_from_name
   implicit none
 
-  integer n, n1, nc, nwf0, m, nwftot, ios
+  integer n, n1, nc, m, nwftot, ios
   logical ok, hf, oep
   !
   !     transform dft in a series of codes for the exchange and
@@ -37,25 +37,10 @@ subroutine ld1_setup
   !
   call set_sl3(sl3,lmx)
   !
-  !  split the spin-orbit states
-  !
-  if (rel == 2) then
-     call occ_spinorb(nwf,nwfx,el,nn,ll,jj,oc,isw)
-     do n=1,nwf
-        oc_old(n)=oc(n)
-     enddo
-  else
-     jj=0.d0
-  endif
-  !
   !  make the correspondence all-electron pseudopotential
-  !  and split the spin-orbit states
   !
   if (iswitch >= 2) then
      do nc=1, nconf
-        if (rel.eq.2) &
-             call occ_spinorb(nwftsc(nc),nwfsx,eltsc(1,nc), &
-             &  nntsc(1,nc),lltsc(1,nc),jjtsc(1,nc),octsc(1,nc),iswtsc(1,nc))
         do n=1,nwftsc(nc)
            nstoaec(n,nc)=0
            do n1=1,nwf
@@ -86,13 +71,9 @@ subroutine ld1_setup
   !
   if (iswitch == 3) then 
      isws=1
-     nwf0=nwfs
-     if (rel==2) &
-          call occ_spinorbps(nwfs,nwfsx,els,nns,lls,jjs,ocs,rcut,rcutus,enls,isws)
      do n=1,nwf
         core_state(n)=.true.
      enddo
-
      do n=1,nwfs
         nstoae(n)=0
         do n1=1,nwf

atomic/ld1inc.f90

@@ -26,7 +26,6 @@ module ld1inc
   real(DP) ::   &
        jj(nwfx),     & ! the total angular momentum
        oc(nwfx),     & ! the occupations of the all-electron atom
-       oc_old(nwfx), & ! saves the occupations of the all-electron 
        zed,          & ! the ionic charge 
        enne,         & ! the number of electrons
        sl3(0:lmx2,0:lmx2,0:lmx2)

atomic/run_test.f90

@@ -16,6 +16,8 @@ subroutine run_test
   use ld1inc
   implicit none
 
+  real(DP) :: oc_old(nwfx)
+
   integer  &
        n, &  ! counter on wavefunctions
        n1,&  ! counter on mesh points
@@ -26,6 +28,10 @@ subroutine run_test
   character(len=1) :: nch
   real(DP) :: dum
 
+  do n=1,nwf
+     oc_old(n)=oc(n)
+  enddo
+
   file_tests = trim(prefix)//'.test'
   open(unit=13, file=file_tests, iostat=ios, err=1111, status='unknown')
 1111 call errore('ld1_setup','opening file_tests',abs(ios))