← Back to branch summary

~nickpapior/siesta/mixing

« back to all changes in this revision

Viewing changes to Src/siesta_analysis.F

  • Committer: Nick Papior
  • Date: 2017-09-08 07:49:32 UTC
  • mfrom: (591.1.41 trunk)
  • Revision ID: nickpapior@gmail.com-20170908074932-l6hiee82h1vjnn2h
Merged trunk-r594-632

Show diffs side-by-side

added added

removed removed

Lines of Context:
Src/siesta_analysis.F
45
45
      use parallel, only: SIESTA_worker
46
46
      use files,       only : label_length
47
47
      use m_energies
48
 
      use m_steps
 
48
      use m_steps, only: final
49
49
      use m_ntm
50
50
      use m_spin,         only: nspin, spinor_dim, h_spin_dim
51
51
      use m_spin,         only: SpOrb, NonCol, SPpol, NoMagn
77
77
      real(dp)          :: polxyz(3,nspin)  ! Autom., small
78
78
      real(dp)          :: polR(3,nspin)    ! Autom., small
79
79
      real(dp)          :: qaux
80
 
      real(dp)          :: Stot             ! Total spin polar.
81
 
      real(dp)          :: Svec(3)          ! Total spin vector
82
80
      real(dp), pointer :: ebk(:,:,:)       ! Band energies
83
81
 
84
82
      integer           :: j, ix, ind, ik, io, ispin
89
87
#endif
90
88
      type(filesOut_t)  :: filesOut  ! blank output file names
91
89
 
92
 
!----------------------------------------------------------------------- BEGIN
 
90
!-----------------------------------------------------------------------BEGIN
93
91
 
94
92
      if (SIESTA_worker) call timer("Analysis",1)
95
93
 
 
94
      ! Check that we are converged in geometry,
 
95
      ! if strictly required,
 
96
      ! before carrying out any analysis
 
97
      
 
98
      quenched_MD =  ( (iquench > 0) .and.
 
99
     $     ((idyn .eq. 1) .or. (idyn .eq. 3)) )
 
100
 
 
101
      if ((nmove .ne. 0) .or. quenched_MD) then
 
102
         if (GeometryMustConverge .and. (.not. relaxd)) then
 
103
            call message("FATAL",
 
104
     $           "GEOM_NOT_CONV: Geometry relaxation not converged")
 
105
            call die("ABNORMAL_TERMINATION")
 
106
         endif
 
107
      endif
 
108
 
96
109
!     All the comments below assume that this compatibility option
97
110
!     is not used.
98
111
!     Note also that full compatibility cannot be guaranteed
157
170
        ! This covers CG and MD-quench (verlet, pr), instances in which 
158
171
        ! "relaxd" is meaningful
159
172
 
160
 
        quenched_MD =  ( (iquench > 0) .and.
161
 
     $                   ((idyn .eq. 1) .or. (idyn .eq. 3)) )
162
 
 
163
173
        if ((nmove .ne. 0) .or. quenched_MD ) then
164
174
 
165
175
          if (relaxd) then  ! xa = xa_last
283
293
         if (IONode) print "(a)", "Writing WFSX for COOP/COHP in " 
284
294
     $                           // trim(wfs_filename)
285
295
         call wwave( no_s, h_spin_dim, spinor_dim, no_u, no_l, maxnh, 
286
 
     .        maxk,
 
296
     .        nkpnt,
287
297
     .        numh, listhptr, listh, H, S, Ef, xijo, indxuo,
288
298
     .        nkpnt, kpoint, no_u, gamma, occtol)
289
299
      endif
367
377
     &     ((isolve.eq.SOLVE_MINIM).and.minim_calc_eigenvalues))
368
378
     &     .and.IOnode)
369
379
     &     call ioeig(eo,ef,neigwanted,nspin,nkpnt,no_u,spinor_dim,
370
 
     &     maxk,kpoint, kweight)
 
380
     &     nkpnt,kpoint, kweight)
371
381
 
372
382
 
373
383
!**   This uses H,S, and xa, as it diagonalizes them again
387
397
!       DMs of the SCF run for the last geometry considered.
388
398
!       This can be considered a feature or a bug.
389
399
 
390
 
!     Print spin polarization
391
 
      if (nspin .ge. 2) then
392
 
        do ispin = 1,nspin
393
 
          qspin(ispin) = 0.0_dp
394
 
          do io = 1,no_l
395
 
            do j = 1,numh(io)
396
 
              ind = listhptr(io)+j
397
 
              qspin(ispin) = qspin(ispin) + Dscf(ind,ispin) * S(ind)
398
 
            enddo
399
 
          enddo
400
 
        enddo
401
 
 
402
 
#ifdef MPI
403
 
!       Global reduction of spin components
404
 
        call globalize_sum(qspin(1:nspin),qtmp(1:nspin))
405
 
        qspin(1:nspin) = qtmp(1:nspin)
406
 
#endif
407
 
        if (nspin .eq. 2) then
408
 
          if (IOnode) then
409
 
            write(6,'(/,a,f12.6)')
410
 
     &       'siesta: Total spin polarization (Qup-Qdown) =', 
411
 
     &       qspin(1) - qspin(2)
412
 
          endif
413
 
          if (cml_p) call cmlAddProperty(xf=mainXML,
414
 
     .         value=qspin(1)-qspin(2), dictref='siesta:qspin',
415
 
     .         units='siestaUnits:spin')
416
 
        elseif (nspin .eq. 4) then
417
 
          call spnvec( nspin, qspin, qaux, Stot, Svec )
418
 
          if (IONode) then
419
 
             write(6,'(a,3f12.6)')
420
 
     &            'siesta: Total spin polarization (x,y,z) = ',
421
 
     &            qspin(3)*2, qspin(4)*2, qspin(1)-qspin(2)
422
 
             write(6,'(a,4f12.6)')
423
 
     &            'siesta: S , {S} = ', Stot, Svec
424
 
            if (cml_p) then
425
 
              call cmlAddProperty(xf=mainXML, value=Stot,
426
 
     .             dictref='siesta:stot',
427
 
     .             units='siestaUnits:spin')
428
 
              call cmlAddProperty(xf=mainXML, value=Svec,
429
 
     .             dictref='siesta:svec',
430
 
     .             units='siestaUnits:spin')
431
 
            endif !cml_p
432
 
          endif
433
 
        endif
434
 
      endif
 
400
      call print_spin(qspin)  ! qspin returned for use below
435
401
 
436
402
!**   This uses the last computed dipole in dhscf during the scf cycle,
437
403
!     which is in fact derived from the "in" DM.