Siesta

Merge lp:~albertog/siesta/4.0-spin into lp:siesta/4.0

  1. 4.0-spin
  2. Merge into rel-4.0
Proposed by Alberto Garcia
Status: Merged
Approved by: Nick Papior
Approved revision: 560
Merged at revision: 562
Proposed branch: lp:~albertog/siesta/4.0-spin
Merge into: lp:siesta/4.0
Diff against target: 502 lines (+195/-88)
15 files modified
Docs/compatibility.tex (+5/-0)
Src/Makefile (+7/-7)
Src/compute_dm.F (+1/-9)
Src/mixer.F (+1/-1)
Src/mulliken.F (+29/-17)
Src/new_dm.F (+1/-43)
Src/print_spin.F90 (+9/-7)
Src/read_options.F90 (+1/-0)
Src/scfconvergence_test.F (+19/-0)
Src/siesta_options.F90 (+2/-1)
Src/state_init.F (+2/-0)
Util/DensityMatrix/README (+4/-0)
Util/DensityMatrix/dm_noncol_sign_flip4.f90 (+108/-0)
Util/DensityMatrix/makefile (+5/-2)
version.info (+1/-1)
To merge this branch: bzr merge lp:~albertog/siesta/4.0-spin
Reviewer Review Type Date Requested Status
Nick Papior Approve
Review via email: mp+342997@code.launchpad.net

Commit message

Spin monitoring in scf cycle. Clarification of sign convention for DM

  * Spin monitoring during scf cycle

  If the fdf variable 'Spin.In.Scf' is set to 'true', the size and
  components of the (total) spin polarization will be printed at every
  scf step. This is analogous to the 'Mulliken.In.Scf' feature.

  The default for the Spin.In.SCF fdf flag is 'true' for calculations
  involving spin.

  A single line of spin information is printed, prefixed by ' spin-moment:'.

  Note that now, for each scf step, the spin and/or mulliken analyses are
  printed *after* the line showing the energies.

  The initial spin polarization for every geometry iteration is also
  printed in 'state_init', replacing the similar calls in 'new_dm'.

  * Use explicitly the same DM_12 sign-convention as in SOC versions

  In the non-collinear case, the code is actually using internally the
  same sign convention for the "up-down" components of H and the DM as
  the spin-orbit-capable (SOC) versions (>= 4.1): The building of the
  dense Hamiltonian has the same form and the off-diagonal components
  of Vxc are identical. Routines diag2g and diag2k are actually using
  the same convention, but with DM* as an intermediate step.

  By changing the interface to the outside world, namely the
  correspondence between spin angles and the DM in new_dm, and the
  calculation of the spin (actually moment) components in spnvec, the
  code can now interoperate (i.e., reuse DMs) with the SOC versions.

  In order to reuse DM information produced by previous versions of
  the program in the non-collinear case, the sign of D(4) of the old
  DM has to be changed. This can be achieved with the
  dm_noncol_sign_flip4 program in Util/DensityMatrix.

To post a comment you must log in.
Revision history for this message
Nick Papior (nickpapior) wrote :

LGTM with the small changes. I think we should just delete the lines that you have commented out (by choosing better alternatives). No need to have them there... :)

review: Approve
Revision history for this message
Alberto Garcia (albertog) wrote :

I agree to remove the old sections.

> On 11 Apr 2018, at 11:41, Nick Papior <email address hidden> wrote:
>
> LGTM

Revision history for this message
Nick Papior (nickpapior) wrote :

Super! Great!

Preview Diff

[H/L] Next/Prev Comment, [J/K] Next/Prev File, [N/P] Next/Prev Hunk
=== modified file 'Docs/compatibility.tex'
--- Docs/compatibility.tex 2018-04-09 13:45:49 +0000
+++ Docs/compatibility.tex 2018-04-11 09:36:55 +0000
@@ -27,6 +27,11 @@
27 internal tables used to compute two-center integrals. Please see27 internal tables used to compute two-center integrals. Please see
28 \opt{Compat.Matel.NRTAB} for details.28 \opt{Compat.Matel.NRTAB} for details.
2929
30 In non-collinear-spin calculations, the sign of the fourth spin
31 component of the density matrix is reversed. To re-use old \code{DM} files,
32 they can be converted using the
33 \code{Util/DensityMatrix/dm_noncol_sign_flip4} program.
34
30 \item[\emph{any} --- 4.0-b2] The following compatibility issues should be remarked when35 \item[\emph{any} --- 4.0-b2] The following compatibility issues should be remarked when
31 comparing with any later version of \siesta.36 comparing with any later version of \siesta.
3237
3338
=== modified file 'Src/Makefile'
--- Src/Makefile 2018-04-09 14:08:42 +0000
+++ Src/Makefile 2018-04-11 09:36:55 +0000
@@ -689,11 +689,10 @@
689mulliken.o: alloc.o atmfuncs.o parallel.o parallelsubs.o precision.o689mulliken.o: alloc.o atmfuncs.o parallel.o parallelsubs.o precision.o
690mulliken.o: siesta_cml.o690mulliken.o: siesta_cml.o
691naefs.o: atmfuncs.o mneighb.o new_matel.o precision.o691naefs.o: atmfuncs.o mneighb.o new_matel.o precision.o
692new_dm.o: alloc.o atomlist.o m_energies.o m_fdf_global.o m_iodm.o m_mpi_utils.o692new_dm.o: alloc.o atomlist.o m_energies.o m_fdf_global.o m_iodm.o m_sparse.o
693new_dm.o: m_sparse.o m_spin.o m_steps.o m_ts_global_vars.o m_ts_iodm.o693new_dm.o: m_spin.o m_steps.o m_ts_global_vars.o m_ts_iodm.o m_ts_options.o
694new_dm.o: m_ts_options.o normalize_dm.o parallel.o parallelsubs.o parsing.o694new_dm.o: normalize_dm.o parallel.o parallelsubs.o parsing.o precision.o
695new_dm.o: precision.o siesta_geom.o siesta_options.o sparse_matrices.o sys.o695new_dm.o: siesta_geom.o siesta_options.o sparse_matrices.o sys.o units.o
696new_dm.o: units.o
697new_matel.o: alloc.o errorf.o interpolation.o matel_registry.o parallel.o696new_matel.o: alloc.o errorf.o interpolation.o matel_registry.o parallel.o
698new_matel.o: precision.o radfft.o siesta_options.o spher_harm.o sys.o697new_matel.o: precision.o radfft.o siesta_options.o spher_harm.o sys.o
699nlefsm.o: alloc.o atmfuncs.o mneighb.o new_matel.o parallel.o parallelsubs.o698nlefsm.o: alloc.o atmfuncs.o mneighb.o new_matel.o parallel.o parallelsubs.o
@@ -768,8 +767,9 @@
768save_density_matrix.o: m_steps.o m_ts_global_vars.o m_ts_iodm.o precision.o767save_density_matrix.o: m_steps.o m_ts_global_vars.o m_ts_iodm.o precision.o
769save_density_matrix.o: siesta_options.o sparse_matrices.o768save_density_matrix.o: siesta_options.o sparse_matrices.o
770savepsi.o: alloc.o parallel.o parallelsubs.o precision.o769savepsi.o: alloc.o parallel.o parallelsubs.o precision.o
771scfconvergence_test.o: m_convergence.o m_energies.o m_wallclock.o parallel.o770scfconvergence_test.o: atomlist.o m_convergence.o m_energies.o m_spin.o
772scfconvergence_test.o: siesta_cml.o siesta_options.o units.o write_subs.o771scfconvergence_test.o: m_wallclock.o parallel.o siesta_cml.o siesta_geom.o
772scfconvergence_test.o: siesta_options.o sparse_matrices.o units.o write_subs.o
773schecomm.o: alloc.o773schecomm.o: alloc.o
774setatomnodes.o: alloc.o parallel.o precision.o spatial.o sys.o774setatomnodes.o: alloc.o parallel.o precision.o spatial.o sys.o
775setspatial.o: alloc.o parallel.o precision.o spatial.o775setspatial.o: alloc.o parallel.o precision.o spatial.o
776776
=== modified file 'Src/compute_dm.F'
--- Src/compute_dm.F 2016-01-25 16:00:16 +0000
+++ Src/compute_dm.F 2018-04-11 09:36:55 +0000
@@ -51,6 +51,7 @@
5151
52 ! e1>e2 to signal that we do not want DOS weights52 ! e1>e2 to signal that we do not want DOS weights
53 real(dp), parameter :: e1 = 1.0_dp, e2 = -1.0_dp53 real(dp), parameter :: e1 = 1.0_dp, e2 = -1.0_dp
54 real(dp) :: dummy_qspin(8)
5455
55!-------------------------------------------------------------- BEGIN56!-------------------------------------------------------------- BEGIN
56 call timer( 'compute_dm', 1 )57 call timer( 'compute_dm', 1 )
@@ -147,15 +148,6 @@
147 endif148 endif
148#endif149#endif
149150
150! Print populations at each SCF step if requested before mixing ......
151 if (muldeb) then
152 if (ionode) write (6,"(/a)")
153 & 'siesta: Mulliken populations (DM_out)'
154 call mulliken( mullipop, nspin, na_u, no_u, maxnh,
155 & numh, listhptr, listh, S, Dscf, isa,
156 & lasto, iaorb, iphorb )
157 endif
158
159! Write orbital indexes. JMS Dec.2009151! Write orbital indexes. JMS Dec.2009
160152
161 if (IOnode .and. iscf==1) then153 if (IOnode .and. iscf==1) then
162154
=== modified file 'Src/mixer.F'
--- Src/mixer.F 2016-01-25 16:00:16 +0000
+++ Src/mixer.F 2018-04-11 09:36:55 +0000
@@ -94,7 +94,7 @@
94 ! entirely correct. It should be moved to the top,94 ! entirely correct. It should be moved to the top,
95 ! or done somewhere else.95 ! or done somewhere else.
9696
97 if (muldeb) then 97 if (muldeb .and. (.not. mixH)) then
98 if (IONode)98 if (IONode)
99 & write (6,"(/a)")99 & write (6,"(/a)")
100 & 'siesta: Mulliken populations after mixing'100 & 'siesta: Mulliken populations after mixing'
101101
=== modified file 'Src/mulliken.F'
--- Src/mulliken.F 2016-01-25 16:00:16 +0000
+++ Src/mulliken.F 2018-04-11 09:36:55 +0000
@@ -533,20 +533,26 @@
533c integer ns : Number of components in spin density matrix533c integer ns : Number of components in spin density matrix
534c real*8 qs(ns) : Spin density matrix elements with the convention534c real*8 qs(ns) : Spin density matrix elements with the convention
535c is=1 => Q11; is=2 => Q22; is=3 => Real(Q12);535c is=1 => Q11; is=2 => Q22; is=3 => Real(Q12);
536c is=4 => Imag(Q12)536c is=4 => -Imag(Q12) (*Note
537c ******* Output *****************************************************537c ******* Output *****************************************************
538c real*8 qt : Total charge538c real*8 qt : Total charge
539c real*8 st : Total spin539c real*8 st : Total spin moment
540c real*8 sv(3) : Spin vector540c real*8 sv(3) : Spin moment vector
541c ********************************************************************541c ********************************************************************
542542!
543! Spin magnetic **moments** in bohr magnetons:
544! m_x = 2 Re{Q12} = 2*qs(3)
545! m_y = -2 Im{Q12} = 2*qs(4)
546! m_z = Q11-Q22 = qs(1)-qs(2)
547!
548! The spin moment for a spin angular momentum of 1/2 (hbar) is (very close to) 1 bohr magneton.
549!
543 use precision550 use precision
544551
545 implicit none552 implicit none
546 integer ns553 integer, intent(in) :: ns
547 real(dp) qs(ns), qt, st, sv(3)554 real(dp), intent(in) :: qs(ns)
548 real(dp) cosph, costh, sinph, sinth, tiny555 real(dp), intent(out) :: qt, st, sv(3)
549 parameter ( tiny = 1.e-12_dp )
550556
551 if (ns .eq. 1) then557 if (ns .eq. 1) then
552 qt = qs(1)558 qt = qs(1)
@@ -561,15 +567,21 @@
561 sv(2) = 0.0_dp567 sv(2) = 0.0_dp
562 sv(3) = st568 sv(3) = st
563 elseif (ns .eq. 4) then569 elseif (ns .eq. 4) then
564 qt = qs(1) + qs(2)570 qt = qs(1) + qs(2)
565 st = sqrt( (qs(1)-qs(2))**2 + 4._dp*(qs(3)**2+qs(4)**2) )571 sv(1) = 2.0_dp * qs(3)
566 costh = ( qs(1) - qs(2) ) / ( st + tiny )572 sv(2) = 2.0_dp * qs(4)
567 sinth = sqrt( 1.0_dp - costh**2 )573 sv(3) = qs(1) - qs(2)
568 cosph = qs(3) / ( sqrt( qs(3)**2 + qs(4)**2 ) + tiny )574 st = sqrt(sv(1)**2+sv(2)**2+sv(3)**2)
569 sinph = -qs(4) / ( sqrt( qs(3)**2 + qs(4)**2 ) + tiny )575!
570 sv(1) = st * sinth * cosph576! ... roundabout way ...
571 sv(2) = st * sinth * sinph577! st = sqrt( (qs(1)-qs(2))**2 + 4._dp*(qs(3)**2+qs(4)**2) )
572 sv(3) = st * costh578! costh = ( qs(1) - qs(2) ) / ( st + tiny )
579! sinth = sqrt( 1.0_dp - costh**2 )
580! cosph = qs(3) / ( sqrt( qs(3)**2 + qs(4)**2 ) + tiny )
581! sinph = qs(4) / ( sqrt( qs(3)**2 + qs(4)**2 ) + tiny )
582! sv(1) = st * sinth * cosph
583! sv(2) = st * sinth * sinph
584! sv(3) = st * costh
573 else585 else
574 write(6,*) 'spnvec: ERROR: invalid argument ns =', ns586 write(6,*) 'spnvec: ERROR: invalid argument ns =', ns
575 return587 return
576588
=== modified file 'Src/new_dm.F'
--- Src/new_dm.F 2016-01-25 16:00:16 +0000
+++ Src/new_dm.F 2018-04-11 09:36:55 +0000
@@ -842,7 +842,7 @@
842 Dscf(ind,1) = (qio + spio * costh) / 2842 Dscf(ind,1) = (qio + spio * costh) / 2
843 Dscf(ind,2) = (qio - spio * costh) / 2843 Dscf(ind,2) = (qio - spio * costh) / 2
844 Dscf(ind,3) = spio * sinth * cosph / 2844 Dscf(ind,3) = spio * sinth * cosph / 2
845 Dscf(ind,4) = - spio * sinth * sinph / 2845 Dscf(ind,4) = spio * sinth * sinph / 2
846 else846 else
847 Dscf(ind,1) = (qio + spio) / 2847 Dscf(ind,1) = (qio + spio) / 2
848 Dscf(ind,2) = (qio - spio) / 2848 Dscf(ind,2) = (qio - spio) / 2
@@ -905,48 +905,6 @@
905905
906 endif906 endif
907907
908 call print_initial_spin()
909
910
911 CONTAINS
912
913 subroutine print_initial_spin()
914 use m_mpi_utils, only: Globalize_sum
915 use sparse_matrices, only: S
916
917 real(dp) :: qspin(nspin)
918 integer :: io, j, ispin, ind
919#ifdef MPI
920 real(dp) :: qtmp(nspin)
921#endif
922
923! Print spin polarization
924 if (nspin .ge. 2) then
925 do ispin = 1,nspin
926 qspin(ispin) = 0.0_dp
927 do io = 1,no_l
928 do j = 1,numh(io)
929 ind = listhptr(io)+j
930 qspin(ispin) = qspin(ispin) + Dscf(ind,ispin) * S(ind)
931 enddo
932 enddo
933 enddo
934
935#ifdef MPI
936! Global reduction of spin components
937 call globalize_sum(qspin(1:nspin),qtmp(1:nspin))
938 qspin(1:nspin) = qtmp(1:nspin)
939#endif
940 if (nspin .eq. 2) then
941 if (IOnode) then
942 write(6,'(/,a,f12.6)')
943 . 'initdm: Initial spin polarization (Qup-Qdown) =',
944 . qspin(1) - qspin(2)
945 endif
946 endif
947 endif
948 end subroutine print_initial_spin
949
950 end subroutine initdm908 end subroutine initdm
951909
952 end module m_new_dm910 end module m_new_dm
953911
=== modified file 'Src/print_spin.F90'
--- Src/print_spin.F90 2017-03-05 16:38:31 +0000
+++ Src/print_spin.F90 2018-04-11 09:36:55 +0000
@@ -45,9 +45,11 @@
45 if (nspin .eq. 2) then45 if (nspin .eq. 2) then
4646
47 if (IOnode) then47 if (IOnode) then
48 write(6,'(/,a,f12.6)') &48 Svec(1) = 0.0_dp
49 'siesta: Total spin polarization (Qup-Qdown) =', &49 Svec(2) = 0.0_dp
50 qspin(1) - qspin(2)50 Svec(3) = qspin(1) - qspin(2)
51 Stot = Svec(3)
52 write(6,'(5x,a,f10.5,2f10.1,f10.5)') 'spin moment: S , {S} = ', Stot, Svec
51 endif53 endif
52 if (cml_p) call cmlAddProperty(xf=mainXML, &54 if (cml_p) call cmlAddProperty(xf=mainXML, &
53 value=qspin(1)-qspin(2), dictref='siesta:stot', &55 value=qspin(1)-qspin(2), dictref='siesta:stot', &
@@ -57,10 +59,10 @@
5759
58 call spnvec( nspin, qspin, qaux, Stot, Svec )60 call spnvec( nspin, qspin, qaux, Stot, Svec )
59 if (IONode) then61 if (IONode) then
60 write(6,'(a,3f12.6)') &62! write(6,'(a,3f12.6)') &
61 'siesta: Total spin polarization (x,y,z) = ', &63! 'siesta: Total spin polarization (x,y,z) = ', &
62 qspin(3)*2, qspin(4)*2, qspin(1)-qspin(2) 64! qspin(3)*2, -qspin(4)*2, qspin(1)-qspin(2)
63 write(6,'(a,4f12.6)') 'siesta: S , {S} = ', Stot, Svec65 write(6,'(5x,a,4f10.5)') 'spin moment: S , {S} = ', Stot, Svec
64 if (cml_p) then66 if (cml_p) then
65 call cmlAddProperty(xf=mainXML, value=Stot, &67 call cmlAddProperty(xf=mainXML, value=Stot, &
66 dictref='siesta:stot', units='siestaUnits:spin')68 dictref='siesta:stot', units='siestaUnits:spin')
6769
=== modified file 'Src/read_options.F90'
--- Src/read_options.F90 2018-04-05 07:11:04 +0000
+++ Src/read_options.F90 2018-04-11 09:36:55 +0000
@@ -1467,6 +1467,7 @@
1467 dm_normalization_tol = fdf_get( 'DM.NormalizationTolerance',1.0d-5)1467 dm_normalization_tol = fdf_get( 'DM.NormalizationTolerance',1.0d-5)
1468 normalize_dm_during_scf= fdf_get( 'DM.NormalizeDuringSCF',.true.)1468 normalize_dm_during_scf= fdf_get( 'DM.NormalizeDuringSCF',.true.)
1469 muldeb = fdf_get( 'MullikenInSCF' , .false.)1469 muldeb = fdf_get( 'MullikenInSCF' , .false.)
1470 spndeb = fdf_get( 'SpinInSCF' , (nspin>1) )
1470 rijmin = fdf_get( 'WarningMinimumAtomicDistance', &1471 rijmin = fdf_get( 'WarningMinimumAtomicDistance', &
1471 1.0_dp, 'Bohr' )1472 1.0_dp, 'Bohr' )
1472 bornz = fdf_get( 'BornCharge' , .false. )1473 bornz = fdf_get( 'BornCharge' , .false. )
14731474
=== modified file 'Src/scfconvergence_test.F'
--- Src/scfconvergence_test.F 2016-01-25 16:00:16 +0000
+++ Src/scfconvergence_test.F 2018-04-11 09:36:55 +0000
@@ -23,6 +23,11 @@
23 use m_convergence, only: converger_t, tolerance23 use m_convergence, only: converger_t, tolerance
24 use m_convergence, only: add_value, is_converged24 use m_convergence, only: add_value, is_converged
2525
26 use m_spin, only: nspin
27 use sparse_matrices, only: Dscf, S, numh, listhptr, listh, maxnh
28 use siesta_geom, only: na_u, isa
29 use atomlist, only: no_u, lasto, iaorb, iphorb
30
26 implicit none31 implicit none
2732
28 integer :: iscf33 integer :: iscf
@@ -32,6 +37,7 @@
32 type(converger_t), intent(inout) :: conv_harris, conv_freeE37 type(converger_t), intent(inout) :: conv_harris, conv_freeE
33 logical, intent(out) :: converged38 logical, intent(out) :: converged
3439
40 real(dp) :: dummy_qspin(8)
35!------------------------------------------------------------------- BEGIN41!------------------------------------------------------------------- BEGIN
3642
37 ! convergence test43 ! convergence test
@@ -59,6 +65,19 @@
59 call wallclock("-------------- end of scf step")65 call wallclock("-------------- end of scf step")
60 endif66 endif
6167
68! Print spin polarization at each SCF step if requested before mixing
69 if (spndeb) then
70 call print_spin(dummy_qspin)
71 endif
72! Print populations at each SCF step if requested before mixing ......
73 if (muldeb) then
74 if (ionode) write (6,"(/a)")
75 & 'siesta: Mulliken populations (DM_out)'
76 call mulliken( mullipop, nspin, na_u, no_u, maxnh,
77 & numh, listhptr, listh, S, Dscf, isa,
78 & lasto, iaorb, iphorb )
79 endif
80
62 converged = .false.81 converged = .false.
6382
64 if (require_energy_convergence) then83 if (require_energy_convergence) then
6584
=== modified file 'Src/siesta_options.F90'
--- Src/siesta_options.F90 2018-04-05 07:11:04 +0000
+++ Src/siesta_options.F90 2018-04-11 09:36:55 +0000
@@ -97,7 +97,8 @@
9797
98 logical :: atmonly ! Set up pseudoatom information only?98 logical :: atmonly ! Set up pseudoatom information only?
99 logical :: harrisfun ! Use Harris functional?99 logical :: harrisfun ! Use Harris functional?
100 logical :: muldeb ! Write Mulliken polpulations at every SCF step?100 logical :: muldeb ! Write Mulliken populations at every SCF step?
101 logical :: spndeb ! Write spin-polarization information at every SCF step?
101 logical :: require_energy_convergence ! free Energy conv. to finish SCF iteration?102 logical :: require_energy_convergence ! free Energy conv. to finish SCF iteration?
102 logical :: require_harris_convergence ! to finish SCF iteration?103 logical :: require_harris_convergence ! to finish SCF iteration?
103 logical :: broyden_optim ! Use Broyden method to optimize geometry?104 logical :: broyden_optim ! Use Broyden method to optimize geometry?
104105
=== modified file 'Src/state_init.F'
--- Src/state_init.F 2018-02-26 15:23:37 +0000
+++ Src/state_init.F 2018-04-11 09:36:55 +0000
@@ -70,6 +70,7 @@
70 logical :: auxchanged ! Auxiliary supercell changed?70 logical :: auxchanged ! Auxiliary supercell changed?
71 logical :: folding, folding171 logical :: folding, folding1
72 logical :: foundxv ! dummy for call to ioxv72 logical :: foundxv ! dummy for call to ioxv
73 real(dp) :: dummy_qspin(8)
73 external :: madelung, timer74 external :: madelung, timer
74 real(dp), external :: volcel75 real(dp), external :: volcel
7576
@@ -332,6 +333,7 @@
332! Initialize density matrix333! Initialize density matrix
333 ! The resizing of Dscf is done inside new_dm334 ! The resizing of Dscf is done inside new_dm
334 call new_dm( auxchanged, numh, listhptr, listh, Dscf)335 call new_dm( auxchanged, numh, listhptr, listh, Dscf)
336 if (nspin > 1) call print_spin(dummy_qspin)
335337
336! Check for size of Pulay auxiliary matrices338! Check for size of Pulay auxiliary matrices
337 call init_pulay_arrays()339 call init_pulay_arrays()
338340
=== modified file 'Util/DensityMatrix/README'
--- Util/DensityMatrix/README 2009-08-20 09:25:37 +0000
+++ Util/DensityMatrix/README 2018-04-11 09:36:55 +0000
@@ -2,3 +2,7 @@
2filter the contribution of certain orbitals to the DM.2filter the contribution of certain orbitals to the DM.
33
4Also, experimental octave scripts to read DM.nc and DMHS.nc files.4Also, experimental octave scripts to read DM.nc and DMHS.nc files.
5
6To allow reuse of non-collinear DM files produced with versions <
74.0.2, use the program dm_noncol_sign_flip4.
8
59
=== added file 'Util/DensityMatrix/dm_noncol_sign_flip4.f90'
--- Util/DensityMatrix/dm_noncol_sign_flip4.f90 1970-01-01 00:00:00 +0000
+++ Util/DensityMatrix/dm_noncol_sign_flip4.f90 2018-04-11 09:36:55 +0000
@@ -0,0 +1,108 @@
1!
2! Copyright (C) 1996-2018 The SIESTA group
3! This file is distributed under the terms of the
4! GNU General Public License: see COPYING in the top directory
5! or http://www.gnu.org/copyleft/gpl.txt.
6! See Docs/Contributors.txt for a list of contributors.
7!
8program dm_noncol_sign_flip4
9
10!
11! Changes the sign of the "4" component of a DM (useful to re-use DM files created
12! with versions prior to 4.0.2 in the non-collinear case).
13!
14! Usage: dm_noncol_sign_flip4
15!
16! The input file must be named "DM" (a symbolic link would do)
17! The output file is called "DMSIGNFLIP4"
18!
19implicit none
20
21integer, parameter :: dp = selected_real_kind(14,100)
22
23integer :: norbs ! Number of atomic orbitals
24integer :: nnzs ! Total number of orbital interactions
25integer :: nspin ! Number of spins
26
27integer :: iog, ispin, iostat, j
28
29integer, dimension(:), allocatable :: numd
30integer, dimension(:), allocatable :: row_pointer
31integer, dimension(:), allocatable :: column
32real(dp), dimension(:,:), allocatable :: dm
33
34!-----------------------------------------------------
35
36open(unit=1,file="DM",form="unformatted",status="old",action="read", &
37 position="rewind",iostat=iostat)
38if (iostat /= 0) then
39 print *, "File DM cannot be opened"
40 STOP
41endif
42
43read(1) norbs, nspin
44print *, "Norbs, nspin: ", norbs, nspin
45
46if (nspin /= 4) then
47 print *, "The DM is not from a non-collinear spin calculation"
48 STOP
49endif
50
51allocate(numd(1:norbs),row_pointer(1:norbs))
52
53read(1) (numd(iog),iog=1,norbs)
54nnzs = sum(numd(1:norbs))
55!
56allocate(column(1:nnzs),dm(1:nnzs,1:nspin))
57
58 ! Compute global row pointer
59 row_pointer(1) = 0
60 do iog=2,norbs
61 row_pointer(iog) = row_pointer(iog-1) + numd(iog-1)
62 enddo
63
64!
65! Column information
66!
67do iog = 1, norbs
68 read(1) (column(row_pointer(iog)+j),j=1,numd(iog))
69enddo
70
71!
72! DM
73!
74 do ispin = 1, nspin
75 do iog = 1, norbs
76 read(1) (dm(row_pointer(iog)+j,ispin),j=1,numd(iog))
77 enddo
78 enddo
79
80 close(1)
81 !
82 ! Change sign
83 !
84 dm(:,4) = -dm(:,4)
85!
86
87open(unit=1,file="DMSIGNFLIP4",form="unformatted",status="new",action="write", &
88 position="rewind",iostat=iostat)
89if (iostat /= 0) then
90 print *, "A file DMSIGNFLIP4 exists. Please remove or change its name"
91 STOP
92endif
93
94write(1) norbs, nspin
95write(1) (numd(iog),iog=1,norbs)
96do iog = 1, norbs
97 write(1) (column(row_pointer(iog)+j),j=1,numd(iog))
98enddo
99do ispin = 1, nspin
100 do iog = 1, norbs
101 write(1) (dm(row_pointer(iog)+j,ispin),j=1,numd(iog))
102 enddo
103enddo
104
105close(1)
106
107
108end program dm_noncol_sign_flip4
0109
=== modified file 'Util/DensityMatrix/makefile'
--- Util/DensityMatrix/makefile 2016-01-25 16:00:16 +0000
+++ Util/DensityMatrix/makefile 2018-04-11 09:36:55 +0000
@@ -21,7 +21,7 @@
21.SUFFIXES: .f .F .o .a .f90 .F9021.SUFFIXES: .f .F .o .a .f90 .F90
22#22#
23default: dm2cdf cdf2dm23default: dm2cdf cdf2dm
24all: dm2cdf cdf2dm dmfilter24all: dm2cdf cdf2dm dmfilter dm_noncol_sign_flip4
25#25#
26OBJDIR=Obj26OBJDIR=Obj
27include ../../$(OBJDIR)/arch.make27include ../../$(OBJDIR)/arch.make
@@ -46,8 +46,11 @@
46dmfilter: dmfilter.o46dmfilter: dmfilter.o
47 $(FC) $(LDFLAGS) -o $@ dmfilter.o47 $(FC) $(LDFLAGS) -o $@ dmfilter.o
48#48#
49dm_noncol_sign_flip4: dm_noncol_sign_flip4.o
50 $(FC) $(LDFLAGS) -o $@ dm_noncol_sign_flip4.o
51#
49clean: 52clean:
50 rm -f *.o dm2cdf cdf2dm dmfilter *.o *.*d53 rm -f *.o dm2cdf cdf2dm dmfilter dm_noncol_sign_flip4 *.o *.*d
51#54#
5255
5356
5457
=== modified file 'version.info'
--- version.info 2018-04-09 14:08:42 +0000
+++ version.info 2018-04-11 09:36:55 +0000
@@ -1,1 +1,1 @@
1siesta-4.0--5611siesta-4.0--561--spin-560

Subscribers

People subscribed via source and target branches