Commit b82a157f authored by giannozz's avatar giannozz

Recently added 'force theorem' calculation of MAE reshuffled:

* generation of noncoliner magnetization from lsda magnetization moved out of
  the routine reading rho into "potinit"
* charge density with noncolinear magnetization so obtained is written to file
  at the end of the noncolinear calculation 
More logical IMHO, no need to tweak low-level I/O routines that can be reverted
to previous state, postprocessing codes read the correct magnetization

Case pw_readfile('rho') deleted: it is useless and confusing (the entire 
routine is a useless wrapper that must disappear, but this another time)


git-svn-id: http://qeforge.qe-forge.org/svn/q-e/trunk/espresso@12170 c92efa57-630b-4861-b058-cf58834340f0
parent 0bc6bc57
......@@ -282,9 +282,7 @@ MODULE io_rho_xml
!
USE io_files, ONLY : tmp_dir, prefix
USE fft_base, ONLY : dfftp
USE spin_orb, ONLY : domag, lforcet
USE noncollin_module, ONLY : angle1, angle2
USE constants, ONLY : PI
USE spin_orb, ONLY : domag
USE io_global, ONLY : ionode
!
IMPLICIT NONE
......@@ -323,30 +321,6 @@ MODULE io_rho_xml
ELSE IF ( nspin == 4 ) THEN
!
IF ( domag ) THEN
!---
! To set up noncollinear m_x,y,z from collinear m_z (AlexS)
!
if(lforcet) then
rho(:,2:4) = 0.D0
write(6,*)
write(6,*) '-----------'
write(6,'(''Spin angles Theta, Phi (degree) = '',2f8.4)') &
angle1(1)/PI*180.d0, angle2(1)/PI*180.d0
write(6,*) '-----------'
file_base = TRIM( dirname ) // '/spin-polarization' // TRIM( ext )
CALL read_rho_xml ( file_base, dfftp%nr1, dfftp%nr2, dfftp%nr3, &
dfftp%nr1x, dfftp%nr2x, dfftp%ipp, dfftp%npp, rho(:,3) )
rho(:,4) = rho(:,3)*cos(angle1(1))
rho(:,2) = rho(:,3)*sin(angle1(1))
rho(:,3) = rho(:,2)*sin(angle2(1))
rho(:,2) = rho(:,2)*cos(angle2(1))
!-------------
else
!
file_base = TRIM( dirname ) // '/magnetization.x' // TRIM( ext )
CALL read_rho_xml ( file_base, dfftp%nr1, dfftp%nr2, dfftp%nr3, &
......@@ -359,7 +333,6 @@ else
file_base = TRIM( dirname ) // '/magnetization.z' // TRIM( ext )
CALL read_rho_xml ( file_base, dfftp%nr1, dfftp%nr2, dfftp%nr3, &
dfftp%nr1x, dfftp%nr2x, dfftp%ipp, dfftp%npp, rho(:,4) )
endif
!
ELSE
!
......
......@@ -894,7 +894,6 @@ interpolate.o : ../../Modules/control_flags.o
interpolate.o : ../../Modules/fft_base.o
interpolate.o : ../../Modules/kind.o
interpolate.o : ../../Modules/recvec.o
io_rho_xml.o : ../../Modules/constants.o
io_rho_xml.o : ../../Modules/fft_base.o
io_rho_xml.o : ../../Modules/funct.o
io_rho_xml.o : ../../Modules/io_files.o
......@@ -1288,7 +1287,6 @@ potinit.o : io_rho_xml.o
potinit.o : ldaU.o
potinit.o : paw_init.o
potinit.o : paw_onecenter.o
potinit.o : pw_restart.o
potinit.o : pwcom.o
potinit.o : scf_mod.o
print_clock_pw.o : ../../Modules/control_flags.o
......@@ -1452,6 +1450,7 @@ read_file.o : ../../Modules/wavefunctions.o
read_file.o : ../../Modules/xml_io_base.o
read_file.o : buffers.o
read_file.o : esm.o
read_file.o : io_rho_xml.o
read_file.o : ldaU.o
read_file.o : newd.o
read_file.o : paw_init.o
......
!
! Copyright (C) 2001-2012 Quantum ESPRESSO group
! Copyright (C) 2001-2016 Quantum ESPRESSO group
! This file is distributed under the terms of the
! GNU General Public License. See the file `License'
! in the root directory of the present distribution,
......@@ -44,11 +44,10 @@ SUBROUTINE potinit()
niter_with_fixed_ns
USE noncollin_module, ONLY : noncolin, report
USE io_files, ONLY : tmp_dir, prefix, input_drho
USE spin_orb, ONLY : domag
USE spin_orb, ONLY : domag, lforcet
USE mp, ONLY : mp_sum
USE mp_bands , ONLY : intra_bgrp_comm
USE io_global, ONLY : ionode, ionode_id
USE pw_restart, ONLY : pw_readfile
USE io_rho_xml, ONLY : read_rho
USE xml_io_base, ONLY : check_file_exst
!
......@@ -86,7 +85,16 @@ SUBROUTINE potinit()
! ... Cases a) and b): the charge density is read from file
! ... this also reads rho%ns if lda+U and rho%bec if PAW
!
CALL pw_readfile( 'rho', ios )
IF ( .NOT.lforcet ) THEN
CALL read_rho ( rho, nspin )
ELSE
!
! ... force theorem: read rho only from lsda calculation,
! ... set noncolinear magnetization from angles
!
CALL read_rho ( rho%of_r, 2 )
CALL nc_magnetization_from_lsda ( dfftp%nnr, nspin, rho%of_r )
END IF
!
IF ( ios /= 0 ) THEN
!
......@@ -126,9 +134,9 @@ SUBROUTINE potinit()
IF (lda_plus_u) THEN
!
IF (noncolin) THEN
CALL init_ns_nc()
CALL init_ns_nc()
ELSE
CALL init_ns()
CALL init_ns()
ENDIF
!
ENDIF
......@@ -255,3 +263,43 @@ SUBROUTINE potinit()
RETURN
!
END SUBROUTINE potinit
!
!-------------
SUBROUTINE nc_magnetization_from_lsda ( nnr, nspin, rho )
!-------------
!
USE kinds, ONLY: dp
USE constants, ONLY: pi
USE io_global, ONLY: stdout
USE noncollin_module, ONLY: angle1, angle2
!
IMPLICIT NONE
INTEGER, INTENT (in):: nnr, nspin
REAL(dp), INTENT (inout):: rho(nnr,nspin)
!---
! set up noncollinear m_x,y,z from collinear m_z (AlexS)
!
WRITE(stdout,*)
WRITE(stdout,*) '-----------'
WRITE(stdout,'("Spin angles Theta, Phi (degree) = ",2f8.4)') &
angle1(1)/PI*180.d0, angle2(1)/PI*180.d0
WRITE(stdout,*) '-----------'
!
! On input, rho(1)=rho_up, rho(2)=rho_down
! Set rho(1)=rho_tot, rho(3)=rho_up-rho_down=magnetization
!
rho(:,3) = rho(:,2)-rho(:,1)
rho(:,1) = rho(:,1)+rho(:,2)
!
! now set rho(2)=magn*sin(theta)*cos(phi) x
! rho(3)=magn*sin(theta)*sin(phi) y
! rho(4)=magn*cos(theta) z
!
rho(:,4) = rho(:,3)*cos(angle1(1))
rho(:,2) = rho(:,3)*sin(angle1(1))
rho(:,3) = rho(:,2)*sin(angle2(1))
rho(:,2) = rho(:,2)*cos(angle2(1))
!
RETURN
!
END SUBROUTINE nc_magnetization_from_lsda
......@@ -284,7 +284,7 @@ MODULE pw_restart
USE ldaU, ONLY : lda_plus_u, lda_plus_u_kind, U_projection, &
Hubbard_lmax, Hubbard_l, Hubbard_U, Hubbard_J, &
Hubbard_alpha, Hubbard_J0, Hubbard_beta
USE spin_orb, ONLY : lspinorb, domag
USE spin_orb, ONLY : lspinorb, domag, lforcet
USE symm_base, ONLY : nrot, nsym, invsym, s, ft, irt, &
t_rev, sname, time_reversal, no_t_rev
USE lsda_mod, ONLY : nspin, isk, lsda, starting_magnetization
......@@ -336,8 +336,10 @@ MODULE pw_restart
CASE( "all" )
!
! ... do not overwrite the scf charge density with a non-scf one
! ... (except in the 'force theorem' calculation of MAE where the
! ... charge density differs from the one read from disk)
!
lrho = lscf
lrho = lscf .OR. lforcet
lwfc = twfcollect
!
CASE( "config" )
......@@ -739,10 +741,9 @@ MODULE pw_restart
! ... CHARGE-DENSITY FILES
!-------------------------------------------------------------------------------
!
! ... do not overwrite the scf charge density with a non-scf one
! ... also writes rho%ns if lda+U and rho%bec if PAW
!
IF ( lscf ) CALL write_rho( rho, nspin )
IF ( lrho ) CALL write_rho( rho, nspin )
!-------------------------------------------------------------------------------
! ... END RESTART SECTIONS
!-------------------------------------------------------------------------------
......@@ -1003,8 +1004,6 @@ MODULE pw_restart
!
lheader = .NOT. qexml_version_init
IF (lheader) need_qexml = .TRUE.
!
ldim = .FALSE.
lcell = .FALSE.
......@@ -1024,8 +1023,6 @@ MODULE pw_restart
lexx = .FALSE.
lesm = .FALSE.
!
SELECT CASE( what )
CASE( 'header' )
!
......@@ -1049,10 +1046,6 @@ MODULE pw_restart
lions = .TRUE.
need_qexml = .TRUE.
!
CASE( 'rho' )
!
lrho = .TRUE.
!
CASE( 'wave' )
!
lpw = .TRUE.
......@@ -1135,6 +1128,10 @@ MODULE pw_restart
lesm = .TRUE.
need_qexml = .TRUE.
!
CASE DEFAULT
!
CALL errore( 'pw_readfile', 'unknown case '//TRIM(what), 1 )
!
END SELECT
!
IF ( .NOT. lheader .AND. .NOT. qexml_version_init) &
......
......@@ -120,6 +120,7 @@ SUBROUTINE read_xml_file_internal(withbs)
USE io_files, ONLY : tmp_dir, prefix, iunpun, nwordwfc, iunwfc
USE noncollin_module, ONLY : noncolin, npol, nspin_lsda, nspin_mag, nspin_gga
USE pw_restart, ONLY : pw_readfile
USE io_rho_xml, ONLY : read_rho
USE read_pseudo_mod, ONLY : readpp
USE xml_io_base, ONLY : pp_check_file
USE uspp, ONLY : becsum
......@@ -301,7 +302,7 @@ SUBROUTINE read_xml_file_internal(withbs)
!
! ... read the charge density
!
CALL pw_readfile( 'rho', ierr )
CALL read_rho( rho, nspin )
!
! ... re-calculate the local part of the pseudopotential vltot
! ... and the core correction charge (if any) - This is done here
......
Markdown is supported
0% or
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment