- force starting charge to be nonnegative

- modified variable-cell dynamics by Cesar da Silva (change if (.false.)
  to if(.true.) in line 823 to enable it)


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

PW/atomic_rho.f90

@@ -48,7 +48,7 @@ subroutine atomic_rho (rhoa, nspina)
   !
   ! local variables
   !
-  real(kind=DP) :: rhoneg, rhorea, rhoima, gx
+  real(kind=DP) :: rhoneg, rhoima, gx
   real(kind=DP), allocatable :: rhocgnt (:), aux (:)
   complex(kind=DP), allocatable :: rhocg (:,:), strf_at(:,:)
   integer :: ir, is, ig, igl, nt, ndm
@@ -149,10 +149,8 @@ subroutine atomic_rho (rhoa, nspina)
      rhoneg = 0.d0
      rhoima = 0.d0
      do ir = 1, nrxx
-        rhorea = DREAL (psic (ir) )
-        rhoneg = rhoneg + min (0.d0, rhorea)
+        rhoneg = rhoneg + MIN (0.d0, DREAL (psic (ir)) )
         rhoima = rhoima + abs (DIMAG (psic (ir) ) )
-        rhoa (ir, is) = rhorea
      enddo
      rhoneg = rhoneg / (nr1 * nr2 * nr3)
      rhoima = rhoima / (nr1 * nr2 * nr3)
@@ -160,9 +158,35 @@ subroutine atomic_rho (rhoa, nspina)
      call reduce (1, rhoneg)
      call reduce (1, rhoima)
 #endif
-     if ( (rhoneg < -1.0d-4) .and. (is == 1 .or. lsda) .or. rhoima > 1.0d-4 ) &
-          WRITE( stdout,'(/"  Warning: negative or imaginary starting charge ",&
-          &2f12.6,i3)') rhoneg, rhoima, is
+     IF ( rhoima > 1.0d-4 ) THEN
+        WRITE( stdout,'(/4x," imaginary charge or magnetization ",&
+          & f12.6," (component ",i1,")  set to zero")') rhoima, is
+     END IF
+     IF ( (is == 1) .OR. lsda ) THEN
+        !
+        !  set starting charge rho or rho up or rho down > 0 
+        !
+        DO ir = 1, nrxx
+           rhoa (ir, is) = MAX (0.d0, DREAL (psic (ir)) )
+        END DO
+        IF ( (rhoneg < -1.0d-4) ) THEN
+           IF ( lsda ) THEN 
+              WRITE( stdout,'(/4x," negative starting charge (spin=",i2,&
+          &                       ") ",f12.6,"  set to zero")') 2*is-3, rhoneg
+           ELSE
+              WRITE( stdout,'(/4x," negative starting charge ", &
+          &          f12.6,"  set to zero")') rhoneg
+           END IF
+        END IF
+     ELSE
+        !
+        !  magnetization for non collinear case (is = 2,3,4)
+        !
+        DO ir = 1, nrxx
+           rhoa (ir, is) = DREAL (psic (ir))
+        END DO
+     END IF
+
   enddo
 
   deallocate (rhocg)

PW/electrons.f90

@@ -709,7 +709,7 @@ SUBROUTINE electrons()
 9018 FORMAT(/'     total magnetization       =',3f9.2,' Bohr mag/cell' &
        &   ,/'     absolute magnetization    =', f9.2,' Bohr mag/cell')
 9020 FORMAT(/'          k =',3F7.4,'     band energies (ev):'/)
-9021 FORMAT(/'          k =',3F7.4,' (',I5,' PWs)   bands (ev):'/)
+9021 FORMAT(/'          k =',3F7.4,' (',I6,' PWs)   bands (ev):'/)
 9030 FORMAT( '  ',8F9.4)
 9032 FORMAT(/'     occupation numbers ')
 9041 FORMAT(/'     the spin up/dw Fermi energies are ',2F10.4,' ev')

PW/vcsubs.f90

@@ -15,6 +15,8 @@ subroutine init (mxdtyp, mxdatm, ntype, natot, rat, ityp, avec, &
   !
   ! rmw (18/8/99)
   !
+  !     Last updated in 04/12/2005 by Cesar RS Silva
+  !
   ! input:
   ! mxdtyp = array dimension for type of atoms
   ! mxdatm = array dimension for atoms (irrespective of type)
@@ -815,12 +817,25 @@ subroutine move (mxdtyp, mxdatm, ntype, ityp, rat, avec, vcell, &
   endif
   if (calc (2:2) .eq.'m') then
      !         WRITE( stdout,109) alpha,nst
+     ! if(.true. ) = original version
+     ! if(.false.) = modified algorithm by SdG
+     !
      if (.false.) then
         do na = 1, natot
            do k = 1, 3
               xx = rat2di (k, na) * rat2d (k, na)
               if (xx.lt.zero) then
                  ratd (k, na) = zero
+
+! ======================================================================
+!        Caution: Under testing!!!!!!!!!
+
+                 rat(k,na)=rat2d(k,na)*rati(k,na)-rat2di(k,na)*rat(k,na)
+                 rat(k,na)=rat(k,na)/(rat2d(k,na)-rat2di(k,na))
+                 rat2d(k,na)=zero
+                 rat2di(k,na)=zero
+! ======================================================================
+
               endif
            enddo
         enddo
@@ -852,8 +867,15 @@ subroutine move (mxdtyp, mxdatm, ntype, ityp, rat, avec, vcell, &
            do l = 1, 3
               xx = avec2d (l, k) * avec2di (l, k)
               if (xx.lt.zero) then
-  !              WRITE( stdout, * ) l, k, avec2d (l, k), avec2di (l, k)
                  avecd (l, k) = zero
+! ======================================================================
+!        Caution: Under testing!!!!!!!!!
+
+                 avec(l, k)=avec2d(l,k)*aveci(l,k)-avec2di(l,k)*avec(l,k)
+                 avec(l, k)=avec(l,k)/(avec2d(l,k)-avec2di(l,k))
+                 avec2d(l,k)=zero
+                 avec2di(l,k)=zero
+! ======================================================================
               endif
            enddo
         enddo